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ABSTRACT 

We model particle growth in a turbulent, viscously evolving protoplanetary nebula, incorporating 
sticking, bouncing, fragmentation, and mass transfer at high speeds. We treat small particles using 
a moments method and large particles using a traditional histogram binning, including a probability 
distribution function of collisional velocities. The fragmentation strength of the particles depends on 
their composition (icy aggregates are stronger than silicate aggregates). The particle opacity, which 
controls the nebula thermal structure, evolves as particles grow and mass redistributes. While growing, 
particles drift radially due to nebula headwind drag. Particles of different compositions evaporate at 
“evaporation fronts” (EFs) where the midplane temperature exceeds their respective evaporation 
temperatures. We track the vapor and solid phases of each component, accounting for advection and 
radial and vertical diffusion. We present characteristic results in evolutions lasting 2 x 10® years. In 
general, (a) mass is transferred from the outer to inner nebula in significant amounts, creating radial 
concentrations of solids at EFs; (b) particle sizes are limited by a combination of fragmentation, 
bouncing, and drift; (c) “lucky” large particles never represent a significant amount of mass; and (d) 
restricted radial zones just outside each EF become compositionally enriched in the associated volatiles. 
We point out implications for mm-submm SEDs and inference of nebula mass, radial banding, the role 
of opacity on new mechanisms for generating turbulence, enrichment of meteorites in heavy oxygen 
isotopes, variable and nonsolar redox conditions, primary accretion of silicate and icy planetesimals, 
and the makeup of Jupiter’s core. 

Keywords: accretion, accretion disks — planets and satellites: formation — protoplanetary disks 


1. INTRODUCTION 


The very early evolution of solids, as they first decouple 
from cosmic gas in the protoplanetary nebula and grow 
into planetesimals such as we see today (asteroids and 
TNOs) remains poorly understood. We call this stage of 
growth primary accretion. The physics of subsequent 
or secondary accretion by mutual collisions (including 
gravitational effects) and sweepup of smaller objects is 
better understood and the process is m ore well charac¬ 
terize d (for a recent review, see, e.g., iRavmond et al.l 
l2014f l. In this stage, gravitational focusing by larger 


objects (so-called runaway accretion) has traditionally 
played a dominant role, and nebula gas a lesser role, but 
recently it has been shown that gas can pl ay an important 
role e ven during this secondary sta ge (|Ormel fc Klahrl 
l201f)l : [l7ambrechts fc Johansenlf2(H^ . Either way, stud¬ 
ies of secondary accretion generally make arbitrary as¬ 
sumptions regarding the sizes of the primary bodies that 

provide their initial condi tions. _ 

Pio ne ering studies by iWeidenschillind (119801 11984 
I1988I 1. iMeakin fc Doniii ( 198811 . Weidenschillin/et al l 


lBhiin&_MunchL (19^), 
iBlum fc WurnJ (2000), 


Donn| 


_ iTml 

Dominik fc Tielend (jlOOTh . . . . . . 

Blum et al.l ( 200^ and others since (see Testi et al. 


2014, and Johansen et al. 2014, 2015 for reviews 
including more recent references) have illustrated the 
importance of the messy physics of stic king and bounc¬ 
ing of grainy, porous, irregular particles. IWeidenschillind 
recognized the important role of gas turbulence 


and the great difficulty of forming large objects by “in¬ 
cremental growth”[3 in the presence of turbulence, due to 
a combination of fragmentation, slow growth, and rapid 
radial drift. Conseq uently, nearly a l l Weidensch i lling’s 
subse quent models (IWeidenschillind I1997L I2000L 12004 
120111) have emphasized conditions where global nebula 
gas turbulence is vanishingly small, allowing a dense 
layer of particles to settle to the midplane; this layer 
itself is s tirred by a small amount of self-generated tur¬ 
bulence (IWeidenschillind 119801: iCuzzi et ai.lll993L [r994 


WeidCTSchilling 1995 : Barrancd 120051 : I Johansen et al.l 
20061 : iBai fc Stonel I2010D . which is however insufficient 


to suppress incremental growth by sticking. 

The possibility of such a nonturbulent environment has 
also motivated a completely different approach to pri¬ 
mary accretion through various kinds of midplane parti- 


^ We define incremental growth as slow and steady growth by 
particles of a given size to larger sizes, by simple physical collision 
and sticking of same-and-smaller-size particles. 
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cle layer instabilities (Safronov 1972, Goldreich & Ward 
1973, Sekiya 1983 et seq., Garaud & Lin 2004, Goodman 
& Pindor 2000, Youdin & Goodman 2005, and subse¬ 
quent workers). In these models, the uncertain prop¬ 
erties of sticking are ultimately made irrelevant by the 
dominant role of local gravitational and/or drag instabil¬ 
ities which lead directly to km-and-larger planetesimals 
(see Ghiang & Youdin 2010 for a review that empha¬ 
sizes work of this type). Midplane instability models all 
require significant initial enhancements of the midplane 
solids-to-gas ratio before they can operate. Because set¬ 
tling of particles into dense midplane layers is frustrated 
by turbulent vertical diffusion, this requirement in turn 
places restrictive conditions on the level of global turbu¬ 
lence and/or global mass enhancement, especially for the 
mm-size particle s that dominate primitive aste roids (for 
a discussion, see iGuzzi fc Weidenschillin^ 1200611 . 

The level of global nebula turbulence remains very 
much a matter of active debate. While the predictions of 
recent MHD models have been trendi ng towards very low 
intensity near-midplane turbulence (iBai fc Stond I20l5 
iBaillMll . purely hyd rodynamical turbu l ence has been 
making a c omeback d^lson et al.l 120131: iMarcns et al.l 


emphasizing the possible differences betw e en yi s cosity 
V and the diffusiyity D. iBirnstiel et al.l (|2010L 1201 il 
l2012al) extended the Brauer results to evolying disks, us¬ 
ing their implicit numerical approach, and added seyeral 
new aspects, including composition-dependent sticking 
(ice vs. silicate); their particle collision yelocities re¬ 
mained single-yalued a s a fu nction of particle size, fol¬ 
lowing [OrmeT^^nii^ (1200711 . Their model opacity did 
not account for the evolying particle size distribution, 
and the only evaporation they included was the silicate 
evaporation front at 1500K where all the solids evapo¬ 
rated. There was no explicit treatment of the behavior 
of vapor of different species. 

In the last decade or so, a large and still growing 
body of experimental work has clarified and refined stick¬ 
ing coefficients for granular particles of various com¬ 
posit ion, porosity, and relative velocity (|Guttler et al.l 
1201011 . The most detailed coagulation models which em¬ 
ploy this detailed sticking physics encounter a “bouncing 
barrier” for mm-cm size particles in moderate turbulence 
(iZsom et al.ll2010ll . This means that, while small parti¬ 
cles can grow by sticking because their collision veloci¬ 
ties are very small, the very act of growing leads them 


for a review). Moreover, a number of meteoritic and 

198C 

; ICarballido et al.l 

2008 

(Mol: IPan & Padoan Mil 

cometary properties point to an extended period of as- 

2013 

. 12014 iHubbardI 

2012 

.1201,31: IPa,n et ail I2ni4a0bll 


teroid formation (iGonnellv et al.ll2012l : iKita fc Ushikubd 
[Ml IKita et all 1201311 and an e xtended radial range 
of mixing of minera l constituents ([Zolenskv et al.l 120061 : 
l.lozwiak et al.llTni2| l. which are both difficult to recon¬ 
cile with a globally nonturbulent nebula in the region of 
planetesimal formation. It seems appropriate to face the 
implications of moderate turbulence squarely, including 
the various barriers it presents to incremental growth. 

Once particles grow into the cm-m size range, de¬ 
pending on location, they drift radially at a rapid pace 
because of headwind drag with the more slowly ro¬ 
tating, pressure-support e d gas (iWeidenschillind 


Cuzzi fc Weidenschiliind 120061 : IWeidenschillind 


1977 


2006 


Johansen et al.l 120141 1. This drift can potentially trans¬ 


port particles great distances while they grow only slowly, 
because turbulence stirs their “feedstock” and leads to 
collision velocities which are erosive or destructive. Be¬ 
cause particles found at any distance from the sun 
may have originated at greater distances, proper mod¬ 
els should ideally treat the entire radial extent of the 
nebula at once. 

Early globa l models with grow th and drift were 
presented by iSteoinski fc ValageasI (| 19961 1199711 and 
iKornet et al.l (1200411. but these w ere simplified in sev- 
eral wavs. iCiesla &: Cuzzi! (1200611 studied the implica¬ 
tions o f drift with simplified growth models, and iGaraudI 
([200711 developed a more comprehensive model includ¬ 
ing an evolving nebula with several volatile speci e s, but 
still with a simplified growth model. iBrauer et al.l (1200811 
were the first to present a complete global (ID) model of 
particle growth and drift with a detailed model of stick¬ 
ing, which was enabled by an elaborate implicit numer¬ 
ical solution to the cumbersome coagulation equation. 
They demonstrated, described, and motivated the vari¬ 
ous barriers to growth (both fragmentation and drift), 
while assuming a s t eady- s tate, non-evolving gas disk. 
iHughes k, Armitag^ (120101 . 1201^ presented models with 
evolving gas disks, but a very simple sticking algorithm 


that exceed the sticking threshold at fairly small sizes 
- in the mm-cm range at 2.5AU. This “bouncing bar¬ 
rier” is in intriguing agreement with the sizes of ubiqui¬ 
tous “chondrules” dominating primitive meteorites, but 
it poses a problem for actually forming the asteroidal 
parents of those very meteorites. Sophisticated vari¬ 
ants of this model showed that weak dm-size agglomer¬ 
ates of dust-rimmed “chondrules” could be formed, but 
only in very low levels of nebula tu rbulence (lOrrn el et ahl 
1200^. IWindmark et al.l (I2012allb[) JGaraud et al.l ([201311 . 
and lPrazkowska et al.l ([2013ll20l4l refined collision mod¬ 
els to include a probability distribution function (PDF) 
for collision velocity, adding a growth path involving frag¬ 
mentation of a sma ller projectile along with mass transfe r 
to a larger target ([Wurm et al.ll200^lBeitz et al.ll2011[l . 
They all found that some very few “lucky” large particles 
could form, although the quantitative details were found 
to be highly sensitive to the numerics of mass histogram 
binning. However, these studies either neglected the ra¬ 
dial headwind drift or assumed arbitrary radial “pressure 
bumps” to prevent it, thus giving only an upper limit on 
the abundance of large “lucky” particles. Our models 
include not only bouncing and partial or total fragmen¬ 
tation effects, mass- and velocity-dependent sticking co¬ 
efficients, velocity PDFs and “lucky” particles, but also 
include full radial drift and a viscously evolving back¬ 
ground gas nebula (section [5]) . 

Sticking properties also depend on particle composi¬ 
tion. From laboratory studies and material properties, 
it is becoming clear that where temperatures are low 
enough for water ice or methanol to form, sticking is sig¬ 
nifica ntly more effective ([Bridges et al.lll99'5lWada et al.l 
so larger and stronger particles can grow more 
robustly (lOknzumi et al.ll2012tl . Thus, accretion might 
operate differently inside and outside of the “snowline”, 
strengthening the need for global models. Growth might 
be faster in the cold outer regions, even while subsequent 
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radial drift might carry this material to the inner solar 
system where the now ice-free refractory material might- 
grow less robustly. Our models include these material- 
dependent sticking properties. 

Of course, the nebula is warmer at smaller radial dis¬ 
tances from the sun, due to a combination of solar heat¬ 
ing and viscous dissipation. Drifting particles do not just 
get “lost t o the sun” but evapo r ate their volatiles along 
the wav (iMorfill fc Volkl 11984 iHueso fc Guillotl 120031: 

Cuzzi fc Zahnlell2004 iKornet et al.ll2004lCiesla fc Cuzzil 

200(11: lOaraudI 1200711 . Each volatile, including ices and 
silicates, has its own Evaporation Front (EF) and our 
model treats them all. Such EFs have a number of impor¬ 
tant implications, creating deviations from some uniform 
“cosmic” abundance of the different materials. Amongst 
the meteoritical implications of such deviations are the 
enrichment of nebula vapor in H^jO, with chemical an d 
mineralogical implications (|Fedkin fc GrossmanI 1200011 . 
and the transport of 0-isotopes formed and frozen out 
in the outer nebula to the inner nebula where they can 
becom e incorporated into meteorites (lYurimoto et al.l 
l2f)0l . Each of these processes has a timescale associ¬ 
ated with it (now poorly known), that can in principle 
be modeled and tied to meteo ritical data. In st udies 
of other mate rials besides HtO. [Pasek et al.l (|2005ll and 
iCieslal (120151) have mode led variable sulfur chemistry, 
and lYang fc Cieslal (|2012ll have modeled temporal and 
radial variations of the nebula D/H ratio due to trans¬ 
port, but in neither case were radial drift of large par¬ 
ticles included. It is a goal of our models to treat these 
processes. 

In calculating the thermal evolution of the nebula, our 
model incorporates two other significant advances. First, 
the nebula opaci ty is c alculated self-consistently as par¬ 
ticles grow (Sec. 12.3.211 . This is very important because 
once a particle exceeds the typical wavelength of thermal 
emission, its opacity decreases linearly with its radius. In 
all other models to date, either ISM (tiny grain) opaci¬ 
ties are assumed even while we know growth must be oc¬ 
curring, or some other arbitrary “representative” value 
of opacity is assumed. Secondly, we model the plausi¬ 
ble luminosity evolution of the early sun. In the first 
10^106 years of its life, the sun’s lumino sity is 10 x to 
3x larger than its main sequence value (|Kusaka et al.l 
II970I : iD’Antona fc Mazzitellil[l994f) . Naturally this will 
have implications for the locations of the “snowline” and 
other evaporation fronts. 

In addition to the bouncing, fragmentation, and ra¬ 
dial drift barriers (which are driven by aerodynamical 
forcing of particle velocities by eddies on a range of 
scales), turbulent nebulae present one additional bar¬ 
rier to incremental growth. In global turbulence, small 
gas density fluctuations associated with pressure fluctu¬ 
ations amongst large scale eddies gravitationally excite 
random velocities of objects in the km- and larger size 
range, much like giant molecular clouds scatter stars in 
our galaxy. These relative velocities actually increase 
with planetesimal size because of the slower damping 
for la r ger object s (iNelson fc Gressel 1201(11: iGressel et al.l 
l20I2ll . 1lda et al.l (j2008[l showed that these velocities are 
strongly disruptive for I — 10 km objects for nominal 
levels of nebula turbulence (objects of 100km size and 
larger are stabilized by their own self gravity). This serial 
gauntlet of barriers to growth in turbulence was discussed 


furthe r bv lOrmel fc Okuzumil (|2013f ) and iJohansen et ahl 

lIMl . 

Recent years have seen the emergence of “leapfrog” 
models that circumvent all these barriers, producing 
lO-lOOkm size primitive bodies directly, in one stage, 


els (ICuzzi et al.ll2001l 

2ooa 

20IC 

,l2014bl ICuzzi & HoffanI 

20121 Johansen et al.l 

20071 

200C 

, 120111: Ichambersll2010l: 


(for a recent review, see l.lohansen et ^120151 ) but they 
share a preference for local enhancement of solids over 
the cosmic abundance ratio by factors of 10 or more. It 
is possible that radial drift can scavenge the outer reaches 
of the nebula, and contribute to such an enhancement of 
solids in the inner nebula; even the region of TNO forma¬ 
tion might, in principle, be augmented in solids by strong 
radial drifts scavenging the r arefied outer nebula at hun¬ 
dreds of AU from the sun (ISteiiinski fc ValageasI 119961 
119971 : iHughes fc Armitag^ I20101 . 1201211 . It may even be 
that such radial redistribution mi ght lead to the “stubby” 
radial distribution of solids that iDesdil (1200711 noted is 
the implication of pre- migration planetary distributions. 
iCiesla fc Cuzzil (|200(ll l found only small degrees of en¬ 
hancement for the inner solar system, but they made 
several simplifying assumptions. In this paper we dis¬ 
cuss the possibilities for large-scale rearrangement of the 
solids-to-gas ratio in the nebula, as a function of time. 

Given the wide range in timescales associated with 
modeling the nebula over its entire radial extent, which 
could extend to as much as ~ 1000 AU, coupled with 
the computationally expensive i mplementation of the 
collisional coagulation equation (iSmoluchowskil 1191611 
has previously made such modeling efforts prohibitive. 
This has typically restricted global models to, for ex¬ 
ample, eithe r studying dust growth and redistr i butio n 
alone (e.g., iBrauer et al.1 120(181 iBirnstiel et al.1 1201011 . 
or studying compositional enhancements at evapora¬ 
tion fronts ( EFs) due to simplified assumptions about 
growth ( e.g.. iCuzzi fc Zahnlil2004 ICiesla fc Cuzzill2006t 
Klaraudl l2007ll or ones that treat growth more car e¬ 
fully, but neglect ra dial drift (jWindmark et al.ll2012al lbl: 
iGaraud et al.ll201^ . Taken together these processes are 
likely to be critical to the evolutionary history of solids 
and condensibles. 

Our model incorporates what may be an ev e n fast er nu¬ 
merical approach than that of IBrauer et al.l (j2008l l (the 
method of moments, see section I2.4|l which avoids direct 
calculations of the detailed size distribution for small par¬ 
ticles, while still continuing to treat the larger particles 
explicitly with detailed size, strength, and relative veloc¬ 
ity distributions, and our opacity uses the local, evolving 
size distribution. It is the large particles that are of most 
interest for radial transport of solids, and it is the details 
of their sizes that may ultimately distinguish between 
various proposed “leapfrog” models of planetesimal for¬ 
mation, so they receive a higher-fidelity treatment using 
the full coagulation equation. 


2. NEBULA MODEL 

The simulations presented in this paper are done us¬ 
ing a new 1 -|- ID radial nebula code that is capable of 
simultaneously treating particle growth and radial migra- 
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Table 1 

Selected Symbols and Parameters used in this paper 


Symbol 

Definition 

Section 

c 

gas sound speed 

sec. 2.1 

^v,d 

(vapor, “dust”) diffusion coefficient 

Sec. 2.2.2 

hd 

effective vertical scale height of “dust” subdisk 

Sec. 2.2.3 

i, j, k 

indices for species, radial bin, and mass 

Sec. 2 

i* 

time variable stellar luminosity 

Sec. 2.3.1 

m, m', rtik 

particle mass 

Sec. 2.4 

m*, r* 

fragmentation mass and radius 

Sec. 2.4.2 

M, 

mass accretion rate, general and for species i vapor and dust 

Sec. 2.1, b.2 

Mu 

initial disk mass 

sec. 2.1 


total mass of “migrators” rtik in a radial bin 

Sec. 2.4.4 

QT: Qd 

Toomre parameters for gas and particle subdisk 

Sec. 2 

Q*, Co 

specific energy for fragmentation and bouncing 

Sec. 2.4.2 

Qa, Qs, Qe, Qe 

radiation transfer efficiencies of particles 

Appendix A.3 

i?, Ro 

radial coordinate, and initial disk radius 

Sec. 1, 2.1 

r, Tk 

particle radii, general or for mass 

Sec. 2.1, 2.4 

ryi 

radius of particle containing most mass in radial bin 

Sec. 2.4.6 


radius of largest particle in radial bin 

Sec. 2.4.1, 3.1.1 

StjStLjStivi 

particle Stokes number = £s^ 

Sec. 2.2.1, 3.2 

St*,Stb,Std 

fragmentation, bouncing and drift Stokes numbers 

Sec. 3.2, 4.2 

5b,*(m, m'); Sk,i 

sticking coefficients for bouncing or fragmentation 

Sec. 2.4.2, 2.4.4 

T, Tpb 

disk midplane and photospheric temperature 

Sec. 2.3.1 

is 

particle stopping time 

Sec. 2.2.1 

li, V 

gas radial and azimuthal velocity 

Sec. 2.4.3 

Uk, Vk 

radial and azimuthal velocity for mass 

Sec. 2.4.3 

yrsA 

radial drift velocity of rrik 

Sec. 2.2.2, 2.4.5, Eq. 51 

Vs 

viscous accretion velocity of nebula gas 

Sec. 2.1, Eqs. 2, 6 

Fv.d 

radial drift velocity (vapor, mass-weighted solid) 

Sec. 2.2.2 


turbulence parameter 

Sec. 2.1 

v.d 

"i 

(vapor, solid) mass fraction in species i 

Sec. 2.2 

^imin? ^isync 

minimum code timestep and synchronization timestep 

Sec. 2 

AVpg 

particle-gas relative velocity 

Sec. 2.4.2, 2.4.3 

AFpp, AVJ 

particle-particle relative velocity 

Sec. 2.4.2, 2.4.3, 2.4.4 

K. 

Rosseland mean opacity 

sec. 2.3.1, 2.3.2, Appendix A.3 


molecular mean free path 

Sec. 2.2.1 


turbulent and molecular viscosity 

Sec. 2.1 

n 

Keplerian orbit frequency 

Sec. 1, 2 

p 

gas mass density 

Sec. 2.1 

Pk 

“dust” particle mass density for mass 

Sec. 2.2.2, 2.4.3 

Pk 

“migrator” particle mass density for mass 

Sec. 2.4.4 

Pp 

particle material (internal) mass density 

Sec. 2.4.1 

Pd 

volume mass density of solids in “dust” population 

Sec. 2.2.2 

cr, cr(m, m') 

collision cross section 

Sec. 2.4.2, 2.4.4 

s 

gas surface mass density 

Sec. 2 


(vapor, solid) surface mass density in species i 

Sec. 2.2.2 

T 

optical depth at thermal wavelengths = acE 

Sec. 2.3.1 


tion over many decades of particle size, while simulating 
the dynamical and thermal evolution of the circumstel- 
lar gas disk. Our model includes self-consistent growth 
and radial drift of particles of all sizes, accounts for the 
vertical diffusion and settling of smaller grains, radial 
diffusion and advection of dust and vapor phases of mul¬ 
tiple species, and a self-consistent calculation of opacity 
and disk temperature which allows us to track the evap¬ 
oration and condensation of the various species as they 
are transported throughout the disk. Our code is paral¬ 
lelized in radial bins which is a natural step in attacking 
problems of this magnitude. The development of further 
innovations and techniques to treat some of the more 
computationally expensive processes helps to make the 
problem even more tractable up to and including addi¬ 
tional parallelization in mass bins. 

In our code, several indices will be used. In general, 
we reserve the index j to refer to quantities that are 
functions of semimajor axis, k will be used to refer to 
the histogram of particle masses, whereas the index i 


will exclusively be used for compositional species. As the 
need arises for additional indices, we assign the index I 
for dummy variables. For our model variables, a lack of 
subscript generally refers to a nebular quantity such as 
the gas surface density S, or the pressure scale height 
H. The subscripts or superscripts d, v or m refer to 
properties of the dust, vapor or migrators as described 
in the following sections. We specifically draw attention 
to the fact that there are many different velocities that 
will be discussed in this paper; the ones most widely used 
are referenced in Table 1. Finally, quantities that can 
be generally referred to will be done so using a subscript, 
such as the total volume density of solid material pd, 
whereas the same quantity that is specific to a radial bin 
will have the subscript transposed to a superscript with 
the index used as the subscript, e.g., p'j. Most of the 
relevant parameters used in our code are summarized in 
Table 1. 

We dehne the minimum time step Atmin in our code 
to be a fraction of the innermost radial bin’s orbital pe- 
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riod P = 2tt/VL^ where 17 = is the orbital 

frequency, M* is the stellar mass, G is the gravitational 
constant, and R is the semi-major axis. This presents a 
dilemma in that the dynamical times in the innermost 
portions of the disk where evolution generally happens 
more quickly can be orders of magnitude shorter than 
those in the outer regions of the disk. Using the same 
time step globally is thus inefficient. 

In order to somewhat circumvent this problem, we em¬ 
ploy an asynchronous time stepping scheme which works 
as follows. Each radial location in the disk R has its 
own time step At associated with it that is the same 
fraction of its orbital period as is Atmin- We compare 
the time “elapsed” in a radial bin with the total sim¬ 
ulation time tsim = I^Atmin, where I is the number of 
iterations the code has currently executed. If a bin has 
executed N steps, then its next execution will occur if 
{N + l)At < /Atinin- When the condition is satisfied, 
the counter N is incremented by 1 and a time step is 
executed for that bin. In this fashion the global “time 
step” for the code is the fraction of the orbital period 
we choose, usually < P/4. The innermost radial bin is 
called every iteration, but we avoid unnecessary calls to 
other radial bins. 

Having each radial location effectively at a different 
evolutionary time poses a problem for the transport of 
material across radial boundaries, either through diffu¬ 
sion or radial drift. We thus introduce the concept of 
the “synchronization time”, a predefined number of steps 
Nsync in which we periodically bring the simulation to 
the same time globally. When MOD(/, Ngync) = 0, all 
radial bins are executed with the time step Atsync = 
AsyncAfinin — NAt. In this paper we choose values 
Nsync A 100. During the synchronization step we exe¬ 
cute g loba l calculations such as the nebula gas evolution 
(Sec. 12.1|) . solve the diffusion-advection equation (Sec. 
IMD, determine the disk temperature (Sec. 12.3.ip and 
write out simulation data. On the other hand, radial 
drift of material can be on time scales which are quite 
fast so that waiting for a synchronization step is not prac¬ 
tical. Instead, we have developed a metho d to a ccount 
for radial drift that is called at all At (Sec. I2.4.5p . Once 
a synchronization step is completed, the time counters 
(/, N) are set back to zero. In full global simulations, 
this synchronization step will still happen much sooner 
than an orbit period in the outermost portions of the 
disk, but the savings in time can be significant. 


2.1. Gas Evolution 

For this work, we use a radial 1-D model for the time- 
dependent evolution of the gas surface density Ii{R,t) 
and radial velocity Vg{R,t). In the case that the grav¬ 
itational potential is due to a central point mass M*, 
the equations for the radial evolution and velocity of the 
nebula gas under viscosity v can be derived from the con¬ 
tinuity and angul ar momentum conservation equations, 
and are given by (|Prin 21 mm) 


dt 


2^ 




( 1 ) 




(2) 


The corresponding disk mass accretion rate is then M = 
—2'KRYjVg where the sign is chosen such that a positive 
value of M indicates accretion onto the central star. The 
total viscosity u can be related to “a—models” in which 
v = atcH = Qftc^/U, where at parametrizes the turbu¬ 
lent intensity, PI is the nebula gas pressure scale height, 
and the gas sound speed can be expressed in terms of 
temperature T as 


c = 



(3) 


In Eq. ([3]) above, k-Q is the Boltzmann constant, the 
adiabatic index 7 = 1.4 for a diatomic gas and /ih = 
3.9 X 10“^^ g is the mean mass per molecule in a mixture 
of hydrogen gas that has ~ 20% helium by number. 

It is generally assumed that the disk is gravitationally 
stable to fragmentation. We monitor the stabili ty of the 
disk i n our code through the Toomre parameter (iToomrd 
[1961 


Qt — 


cU 

ttGS' 


(4) 


Values of Qt less than unity lead to disk fragmentation, 
but for values Qt ^ 2 the disk is described as weakly 
unstable and may lead to the formation of clumps. How¬ 
ever, these c lumps only fo rm if the disk is atypically mas¬ 
sive or cold (lR afikovl2005fl , and unless there are processes 
to keep the disk unstable, it is assumed that weak grav¬ 
itational instabilities quickly lead to stabilization of the 
disk via the excitation of spiral density waves. The pro¬ 
cess is thus self-limiting to the extent that these waves 
carry away angular momentum that spread the disk, low¬ 
ering E. If Toomre-unstable conditions arise, we can in¬ 
troduce a moderately large value of at until stable condi¬ 
tions resume. In practice, however, this does not happen 
in the disk models presented in this paper. 

We derive initial conditions fo r our disk models using 
the an alytical expressio ns from iLvnden-B e ll fc Pringle! 
dHH), as generalized bv iHartmann et al.l (Il998il . which 
are parametrized by some initial disk mass Md and ra¬ 
dius Rq. In determining the initial gas surface density 
and radial velocity, the value and radial dependence of 
the turbulent viscosity is expressed as a general power 
law o f the form n = vq{R/RqY (e-g-) iHartmann et al.l 
1199811 . These assumptions lead to simple 1-D, vertically 
integrated and a veraged expressions whi ch can be read¬ 
ily derived from IHartmann et al.l (1199811 given here for 
t = 0: 


E(i?,0) 


Md [2- 


ttRI 


R 

Rq 


AL j p-iR/Rof 


(5) 


-li (I 




2-/5 


( 6 ) 

where vq = atCo/Dg, with all quantities evaluated at Rq. 
Typical ranges for the power law exponent that charac¬ 
terize plausible ext remes of nebula r a dial variation are 
( 3 < /3 < 3/2 (e.s.. iCuzzi et al.ll2003ll . IHartmann et al.l 
([19^ favor Md = 0.2 M©, Rg = 10 AU and /3 = 1 
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based on young star statistics and we use these as fidu¬ 
cial values in this paper. The scale parameter Rq can 
also be cho sen to match sola r system specific angular 
mom entum (iCuz zi et aLll ^OO.Il Rq =4.5) or by other cri¬ 
teria (|Yang fc Cieslall2012L Rn =53). Indeed the nominal 
value of Rq has more often been set at larger values 20- 
60AU (iCiesla fc Cuzzill2006l iGaraudI 120071: iBrauer et al.l 
120081 : iHughes fc Armitaed l2012t lYang fc Cieslal l2012il 
than chosen here but the rationale is not always clear. 
We note that these choices simply provide initial condi¬ 
tions for the disk’s surface density profile in terms of the 
initial disk mass Mu, and that the viscosity of the disk 
will in general not follow a single power law distribution 
- either initially, or as a function of time due to particle 
growth, changing opacity, and our sel f-con sistent calcula¬ 
tion of the disk temperature (see Sec. 12.31) . Thus, at sub¬ 
sequent timesteps the actual viscous evolution equations 
are integrated using the actual local properties. How¬ 
ever, as noted in section 01 these choices do influence 
disk evolution at early times of interest. 

2.2. Evolution of Dust and Vapor 

A self consistent calculation of the disk temperature 
requires that we also know the distribution of dust and 
vapor within the disk, and the initial distribution is es¬ 
tablished during our first calculation of the disk temper¬ 
ature. Our model is capable of tracking any number of 
species in both the vapor and solid phase. We define the 
concentrations of each species i such that 

S v,d 


is the fractional mass of constituent i, or the ratio of 
surface density of solids or vapor to the gas surface den¬ 
sity. In this paper, we include only five different species 
which are listed in Table 2 along with their correspond¬ 
ing condensation temperatures R, densities pi and ini¬ 
tial concentrations Xi in the solid state. We denote by 
fd, fv and /in the total fractional amounts of the dust, 
vapor, and migrator (see Sec. 12.4.41) components within 
the disk. For instance, /d = These quantities are 

updated locally at every time step At, and globally at 
every synchronization step tsync- 

2 . 2 . 1 . Particle Stokes Number and Stopping Times 

Treating the growth and redistribution of solids in the 
nebula concomitant with its gas and thermal evolution 
is a key aspect of this work, and is essential for a self- 
consistent model of the nebula. Initially, the dust-to¬ 
gas ratio in the disk is assumed to be of order ~ 10“^ 
which is consistent with values for the ISM. At this av¬ 
erage abundance, the mass volume density of solid ma¬ 
terial pd does not affect the motion of the gas unless 
settling of dust particles leads to a layer in the midplane 
with density Pd > p, where p is the gas mass density 
(|Nakagawa et al.iri986f) . In fact, we can define an equiv¬ 
alent condition to Eq. o for the dust Qd = Df /irCpd, 
such that if Qd ^ Qt, the condition for gra vitational 
instability of the dust layer may be satisfied (|Safronovl 
119911: ICuzzi et"^ 1199311 . However, even for Qd ^ Qt, 
collectively the dust can play a significant role in the 
disk’s thermal evolution and thus affect E and Vg through 
the temperature distribution, if a large fraction of the 


dust particle sizes remain small such that the opacity re¬ 
mains high (Sec. 12.3.21) . In particular, differences in the 
growth and migration rate, as well as the possibility for 
evaporation and condensation of solid grains, can lead to 
significant radial variation in the distribution of solids. 

For a large range of solids-to-gas ratios, the influence of 
the nebula gas on the motion of the particles over the full 
spectrum of sizes can be described by the dimensionless 
Stokes number 


St ^ (8) 

^ed 

where tg is the particle stopping time, and ted is some 
eddy turnover time, which we choose here to be the in¬ 
tegral scale, or turnover time of the largest eddy in the 
turbulence. This is generally assumed (and has been 
shown) to be for global turbulenc e (ICuzz i et al.l 

120011 : 1.Tohansen et al.l[2007l : iCarballido et al.il2010() . The 
particle stopping time is defined as the time needed for 
the gas drag force to dissipate a particle’s momentum rel¬ 
ative to the gas. The drag force on a particle of radius r 
depends on the size of the particle relative to the molec¬ 
ular mean free path Amfp, and can be separated into two 
distinct flow regimes: 


ts = —- if ?■< (9/4)Amfp; (9) 

P ^ 


ts 


8 Pp r 
3 p CdAFpg 


if r > (9/4)Amfp, 


( 10 ) 


where Pp is the particle material density fsection |2.4.1|) . 
Equation [9] describes the Epstein flow regime in which 
smaller grains are well coupled to the gas flow and rel¬ 
ative velocities between grains are small. In the larger 
particle Stokes flow regime (Eq. [10]), particles become 
increasingly less coupled to the gas flow and their stop¬ 
ping times are affected by the drag coefficient Cd which 
depends on a particle Reynolds number Rep that is it¬ 
self a function of the relative velocity AWp- between the 
particle and the gas ()Weidenschilling|[l977lf . 

In our code, we calculate the particle-to-gas and 
particle-to-particle relative velocities over the evolving 
size distribution which can cover many decades in mass. 
As growth proceeds to larger sizes, some particles will re¬ 
main in the Epstein flow regime, while larger ones will be 
subject to Stokes flow. We use a bridging expression that 
pr ovides a smooth transition between the two regimes af¬ 
ter [PodoiaketXn (1198811 . Our treatment of the stopping 
times is described in further detail in Appendix A.4. 


2 . 2 . 2 . Radial Diffusion-Advection 

We determine the radial motion of the solid and va¬ 
por fractions of all compositional species by solving the 
advection-diffusion equation 


dt 


RdR 


RUv.dS- 


dR 


■ — -n'Vv.d^ 


( 11 ) 

where is the surface density for dust (d) or vapor 
(v) of species i, lA.d is a net, inertial space advection 
velocity, and D^ d is the diffusivity. 
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Table 2 

List of Species 


Species 

0. 

bf 

psp (g cm 3) 

^sp (T < Tsp) 

Iron 

1810 

7.8 

1.26 X 10-* 

Silicates 

1450 

3.4 

3.41 X 10-3 

Troilite 

680 

4.8 

7.68 X 10-4 

Organics 

425 

1.5 

4.132 X 10-3 

Waterice 

160 

0.9 

5.55 X 10-3 


The separately treated term Si represents sources and 
sinks for “dust” and vapor for species i which, in our 
treatment, includes the growth, radial drift and destruc¬ 
tion of migrating material (Sec. I2.4.4|l . The sign of the 
advection velocity in Eq. m is such that 14 ,d < 0 indi¬ 
cates inward radial velocity, and 14,d > 0 outward radial 
velocity. 

For the vapor phase, w e assume Dy = Dg = v/Scg 
((Hughes fc Armitagell2010ll . where the gas diffusivity Dg 
is the ratio of the viscosity to the Schmidt number Scg 
(which we take to be unity in this work, see Appendix 
B) and the vapor advection velocity is just that of the 
gas; Vy = Vg (equations [D and (5]). The particle dif- 
fusivity Uh = .0^/(1 -I-St^) ((Youdin fc Lithwi^ 120071 : 
iCarbailido et al.l[20ilh and 14 are determined at radius 
R through a mass-weighted mean of all grain sizes in the 
dust population (indexed by k): 



f V 

(12) 

Vi + st^J ’ 

II 

Pk T/rad 

- 5 

Pd 

(13) 


where pd = Pk total mass volum e density of 

solids in the subdisk layer (see Sec. 12.2.31) . The radial 
velocity in Eq. dEl) properly takes into account 

the radial drift of particles rela tive to the gas due to 
the local pressure gr adient le.g.. iNakagawa et aklllQSfl 
iTakeuchi fc Linll2002ll and is described in more detail in 
Sec. 12.4.^ (Eg. [50(1. 

The sign of the mean radial drift velocity depends on 
the particle size (large particles drift inwards and small 
ones are advected outwards). We utilize a power law 
distribution in particle mass for the dust population with 
a typical exponent q = 11/6, which is representative of a 
collisional population (see Sec. I2.4.1|) . so that most of the 
mass will be in the largest sizes and the mean velocity will 
always be inward. However, the mass average approach 
by itself as written in Eq. (m would fail to take into 
account that the smallest grains would still move with 
the gas, and in the outer parts of the disk, this flow 
may be outward. In order to account for this in our 
code, we treat the dust as two separate populations and 
define a mean velocity for each, which characterizes 
the mass fraction in small particles that move with the 
gas, and V/" which characterizes the larger grains that 
radially drift inward. Before calculating the mean drift 
velocities, we determine the grain size that separates the 
two populations, and define the total fractional mass in 
each population which is used as a weighting factor when 


solving the advection-diffusion equation (see Appendix 
A.2). 

2.2.3. Vertical Diffusion and Subdisk Height 

We treat the vertical diffusion and settling of dust 
grains in the “-|-1D” part of our code using an analyt i¬ 
cal solution combining elements of iDubrulle et al.l (1199511: 
iCuzzi fc Weidenschillingl ((200611 and lYoudin fc LithwickI 
((2007(1 to calculate the particle distribution as a function 
of height z in the disk. The vertical scale height of 
particles with mass m is defined as 

hdim) = + (14) 

This solution bears a resemblance to that of IGaraudI 
(I2007L her Eq. 22) but includes a number of subtleties 
(see Appendix B). 

To calculate a characteristic scale height for the 
local particle subdisk as a whole, we specify a represen¬ 
tative particle mass as either half the fragmentation bar¬ 
rier mass TO, or the largest particle mass tol/2 (if no 
particles have yet reached to*; see Sec. I2.4.2L based on 
focused 2D calculations that follow growth as a function 
of height at a given radius R. In future refinements, a 
value of /id could be tracked for all particle sizes, giving a 
complete size distribution as a function of altitude (as for 
instance in Sec. 12.4.31 Eq. [38]). The relative velocities 
(Sec. I2.4.3P are also defined at the same representative 
/id- The layer of thickness /id defines the volume acces¬ 
sible to particles growing in the midplane. 

2.3. Disk Thermal Evolution 
2.3.1. Temperature 

We calculate the protoplanetary disk midplane and 
photosphere temperatures self-consistently, assuming 
that the nebula is heated by a combination of internal 
viscous dissipation E^, at a rate proportional to E and 
V and thus primarily near the midplane, and external il¬ 
lumination A* by the stellar luminosity L* (which can 
vary with time; see below). The stellar luminosity heats 
the upper layers of the disk on b oth sides, and thus indi¬ 
rectly the material beneath (e.g..lRuden fc Pollackll 19911 : 
iWoolum fc Cassenlll999t iCieslall^lOH . The thermal en¬ 
ergy of the disk is radiated into space from the photo¬ 
sphere at Tph{R,t), the altitude where the optical depth 
at thermal wavelengths measured vertically outwards is 
roughly unity. For cosmic abundance and a standard 
ISM-MRN grain size distribution, the photosphere lies at 
a distance iJph ^ ± 3 — 5 H from the midplane. On each 
face of the nebula, we can express the energy balance as 
dsBTpj, = E^/2 + E-i,, where ctsb is the Stefan-Boltzmann 
constant. I n reality, the disk has an optically thin hot 

S here ((Chiang fc Goldreichl Il997t iDullemond et al.l 
which indirectly warms the disk photosphere to 
Tph, but it can be shown that modeling direct deposition 
of solar energy at, and thermal radiation by, the disk 
p hotosphere leads to the same Tph. Note that Eq. (12a) 
of iChiang fc Goldreichl (11997( 1 has an extraneous factor 
of 2^/“^ which is removed by integration over all angles 
of the optically thin emission fr om the superheated ex¬ 
osphe re slab, as for instance in iNakamoto fc Nakagawal 

(HSll). 
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We are primarily interested in the midplane tempera¬ 
ture T, because planetesimals and boulder-size drifting 
particles which transport solids radially and feed EFs 
lie mostly near the midplane and not at high altitudes. 
The midplane temperature is influenced both by and 
E*. The external irradiation is both deposited and re¬ 
radiated at the disk photosphere, producing zero net ver¬ 
tical flux through the rest of the disk in steady state and 
thus a vertically constant temperature below the pho¬ 
tosphere in the absence of other energy sources. The 
energy produced by viscous dissipation must flow verti¬ 
cally away from the midplane (assuming it is primarily 
produced there, as we do) to the photosphere to be radi¬ 
ated away, leading to a vertical thermal gradient. The re¬ 
sulting temperature distribution depends on whether the 
disk i s optically thick or thin. iNakamoto fc Nakagawal 
(jl994l l suggest a “bridging” expression that covers both 
of these regimes (their Eq. A15), which can be simplified 
to: 


asBT^- + (15) 

where r is the full optical depth of the nebula at thermal 
wavelengths. Between ±i/ph, t = kE/2, where k is the 
average thermal opacity. Energy from in f all sh ocks was 
also included by INakamoto fc Nakagawal (119941 1 , but we 
do not treat it in this paper. 

Comparing Eq. (USD to that for the disk surface shows 
that the midplane temperature responds differently to 
internal and external energy sources. Regarding inter¬ 
nal sources, most preyious workers have treated only the 
optically thic k lim it (3r/8 term above). F or instance, 
lCassenl(jl994l l. and lWoolum fc CassenI(Il999f l generalized 
the classical radiative transfer solution for a layer of op¬ 
tical depth r ^ 1 which they defined from the midplane 
outwards as = (3/4)rTj£. This is the classical Ed¬ 
dington solution where all the energy is produced at the 
bottom of the layer, and is thus only half of the value in 
Eq. (fTCl) . Cassen (1993; and others subsequently) merely 
stated that the leading factor of 3/8 applies to the more 
realistic situation where the energy production is verti¬ 
cally distributed and roughly proportional to the mass 
density. The derivation of the 3/8 factor for the regime 
of int erest can be extracted from iShakura fc SunvaevI 
and assumes constant k which, for large optical 
depths, is the Rosseland mean opacity (see Sec. 12.3.21) . 

The 1/2t term in Eq. (TTSl) treats the opposite limit 
where significant internally generated energy must be 
radiated away, but the nebula opacity capable of doing 
this is low, for instance, if most of the grains have evap¬ 
orated. In this optically thin regime the emitted flux 
from each face is given by 2rcrsB’r'^, where r is the full 
optical depth, and the factor of 2 comes from integrating 
intensity over all angles. In this regime, temperature is 
roughly independent of altitude (T ^ Tph). 

We assume a local, vertically integrated viscous diss- 
pation rate Ei, = (9/4)S]zzn^ which is closely relate d 
to the mass accret i on ra te (|Lin fc Papaloizoul II985I1 . 
IShakura fc Sunva"^ (|I973ll note that the local energy 
production rate E^, differs by a factor of three from the 
local rate of release of gravitational potential energy. The 
stellar flux on each disk face is given by 


i? j.- rri f 


(16) 


where is some grazing incidence angle depending 
on disk geometry, and the stellar luminosity is L* = 
AnR^asBT^, with i?* and T* the stellar radius and photo- 
spheric temperature. The incidence angle which deter¬ 
mines the solar disk heating differs considerably between 
so-called “flat” disks, having photospheric height Eph a 
constant fraction of the distance R from the star, and 
“flared” disks where dHpi^/dR increases with i?, which 
leads to much larger thickness. Whether the disk is 
flared or flat has implications for the gas mass densi¬ 
ties we use to determine particle interactions, growth 
and drift (Sec. I2.4.ip . For the more realistic flared 
disks, value s of di were derived by iKenvon fc Hartmaiml 
(ITMl. and iRuden fc PollackI (Il99lh . and reiterated bv 
iChiang fc Goldreichl ( I997I1 . who also derive the general 
radial variation 


R-k ^ d f DIph 


^(^).0.4^+R- ^ 


(17) 


0.005RXu + 0.057?^/’' 


•■AU) 


where Rav is radial distance in AU, and the first term 
on th e RHS is the “flat disk” value le.g.. lAdams fc Shul 
II986I1 . The second term on the RHS can be understood 
as a manipulated version of the incidence angle at the 
photosphere of a flared disk given more obviously by 
dH/dR— E[/R. Thus the illumination term is dominated 
by flared disk geometry in general. Iterative treatment 
of 4>{R), depending on Tph a nd L^, is left for futur e work . 

I n one subtle difference, [ Chiang fc Goldrei^ (119971) 
and IRuden fc PollackI (jI99If) appear to assume th at each 
disk face se es only half of the stellar flux, while iCieslal 
assumes that the entire star is visible and 
the flux normal to itself at the disk face is simply 
L*/47ri?^ = asBT^{R*/R)^- For the entire star to be vis¬ 
ible from some point at (i?, i7ph), the opaque disk must 
vanish inside of a radius Rmin{R), where the photosphere 
height is Eph, min, such that 


77ph, min 4” Rk Rph 4” Rk 


Rn 


R 


(18) 


so using Hpi^/R ^ 0.27?^^ (jChiang fc Goldreichl Il997l l 
at both R and 7?min, we find 


7?min(7?) = R 


D 

I-f 0.2—77 


Rk 


2/7 

i-AU 


1 -1 


(19) 


That is, even from a point on the disk photosphere at 
3 AU from the star, the entire stellar disk can be seen 
unless the opaque nebula disk extends further inwards 
than 37?*, and the stellar disk becomes more visible at 
larger distances. Thus we will assume the full stellar flux 
illuminates each face of the disk (this is possible because 
the disk is flared and not planar). On the other hand, the 
value of (j) adopted bv iCieslal (|2009l[20I0ll is several times 
smaller than the values given by Eq. (El) above, reduc¬ 
ing the stellar flux accordingly. The equation determin- 
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ing the midplane temperature must be solved iteratively 
using a root finding algorithm; 


asBT^ = 

O 





( 20 ) 


where the Rosseland mean opacity k is a function of 
T through the t— dependent (T-dependent) dominant 
wavelength, the evolving particle size distribution, and 
the solids fractions of all species i. 

We em ploy a time variable luminosit y R*(t) using the 
model of lD’Antona &: Mazzitellil (|1994 Table 3) for a 1 
Mq star. We fit a polynomial to these authors’ tabulated 
yalues, which coyer a time scale of 7 x 10“^ — 1 x 10® years 
after collapse. This giyes an initial stellar luminosity of 
L* « 12 L© at 7 X 10^ years, which we use as the starting 
time for our simulations, dropping to perhaps 3 L© at 10° 
years. 


2.3.2. Rosseland Mean Opacity 

The expressions from Sec. 12.3.11 relating the mid¬ 
plane temperature to the sources of energy assume a 
wayelength-independent or grey opacity. The transfer 
of thermal radiation in regions of high optical depth is 
maximized at wavelengths where opacity is low. The 
standard treatment is to define a Rosseland mean opac¬ 
ity K from the basic wavelength-dependent opacity k\, 
weighting the inverse (the transparency) at wavelength 
A by the derivative of the Planck function dB\{T)/dT 
which manifests the local flux gradient: 


1 ^ I ^ TT 

f^,dX dasBT® 


dT 


dX. 


( 21 ) 


The Planck opacity, which is preferred for low o ptical 
depth regions fe.g.. iNakamoto fc Nakaga^ 1 199411 . can 
be calculated similarly through a straight average of the 
wavelength-dependent o pacities kx as weigh ted by the 
Planck function Bx{T). iPollack et al.l (Il994f) show that 
the distinction between the Rosseland and Planck opac¬ 
ities is not large for small solid particles, and thus as 
previously stated we include only the Rosseland mean 
opacity in our temperature calculations. 

We utilize a new opacity model in order to determine 
the K\, which is fully described in Cuzzi et al. (2014; 
their Appendix A contains a derivation of Eq. [H]), 
which can be easily incorporated into evolutionary mod¬ 
els at little computational cost. We utilize realistic mate¬ 
rial refractive indices for a cosmic abundance suite that 
likely characterizes nebula solids: water ice, silicates, re¬ 
fractory organics, iron sulfide and metallic iron (sum¬ 
marized in Table 2). These indices a nd re lative abun¬ 
dances are taken from IPollack et aP (|1994D . but alter- 
nate tabulations ca n readily be used in our code (e.g., 
iDraine fc Led 11984 iHenning et al.l 199911. We briefly 
summarize how the opacity model of iCuzzi et al.l (|2014ari 
is applied in our model in Appendix A.3. 


2.4. Solid Body Growth 

Global modeling of the aggregation and radial evolu¬ 
tion of solids in the protoplanetary nebula is required 
over time scales of millions of years in order to un¬ 
derstand many key aspects of primitive bodies. A key 


feature of our model is the capability of modeling par¬ 
ticle growth over many decades of particle size, from 
sub-micron-sized dust to meter-sized and larger boulders 
which can radially drift large distances as they grow. 

Growth by sticking in the protoplanetary disk starts 
with sub-micron grains which are dynamically coupled 
to the nebula gas, and proceeds incrementally through 
larger sizes that collide at larger relative velocities (Sec. 
12.4.31) . The growth of aggregates continues at a rate de¬ 
termined by local nebula conditions until some “bounc¬ 
ing barrier” is reached where collisional energy cannot be 
dissipated, but remains inadequate to destroy the collid¬ 
ing particles (see Sec. [T]). The particle size where this 
occurs depends on assumed particle strength and the lo¬ 
cal value of turbulent at ■ 

The bouncing barrier is not impermeable, but merely 
slows growth by restricting collision partners. Growth 
beyond the bouncing barrier may continue by accre¬ 
tion of sufficiently smaller particles, through a transi¬ 
tion regime of “sub-migrators” where particles are highly 
susceptible to mutual collisional destruction (fragmenta¬ 
tion) as they drift radially, to “migrators” which have 
grown large enough to have a much lower probability of 
destruction (|Estrada fc Guzzill2009[ l. The latter particles 
can drift large distances and grow further, perhaps even 
into planetesimals. 

The treatment of growth up to the bouncing and/or 
fragmentation barriers has until recently presented the 
most significant challenge, because a full-scale solution 
to the problem of dust coagulation for a size spectrum 
at every spatial location in the disk, and as a function 
of time, has been computationally prohibitive, and as 
a result previous models have been limited in different 
ways (see Sec. [T]). Growth through the transitional sub¬ 
migrator and migrator regimes between dust and plan¬ 
etesimals, and associated radial drift over long times, is 
even less thoroughly studied. In our model, we use the 
moments method of lEstrada fc Guzzil (I2008D to greatly 
accelerate coagulation modeling up to t he bouncing bar¬ 
rier, as described below in Sec. 12.4.11 Growth beyond 
the bouncing barrier through the transitional regimes is 
treated explicitly in the traditional way, as described in 

Sec. mu 

In this paper we do not treat collective concen¬ 
tration or sweepup effects such as streaming insta- 
bilities (|Goodman fc Pindoil 120001 : lYoudin fc Goodni^ 

l2?M see howeve r, _ Sec 14.2)1 . turbulent concentra- 

tion dGuzzi et al.l 120081 1201011. or “pebble accretioii ” 
(|Ormel fc Klahfi 120101 : iLambrechts fc JohansenI l2012fl . 
All these subsequent processes depend strongly on initial 
local conditions (primarily, particle size and abundance) 
which are determined by sticking for some given at in 
the presence of drift, and it is those conditions which are 
the outcome of the models described here. 


2.4.1. Coagulative Grain Growth 

The standard approach to modeling dust coagulation 
involves s olving some form of the collisional coagulation 
equation (iSmoluchowskilfTgT^ 
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df{m,t) _ 1 
dt 2 


pm 

/ K{rn — m',rn')f{rn — 7n\t)x 

Jo 

poc 

f {m',t) dm' — / K{m,m')f{m,t)f{m',t)dm' 

Jo 


, ( 22 ) 


Under the assumption of a powerlaw with fixed expo¬ 
nent q, we can derive an equation for the time rate of 
change of the particle mass mn that characterizes the 
largest mass in the distribution, until tol gets as large 
as the fragmentation ma ss (see next section) (Eq. (28) 
of lEstrada fc Cuz'^l2008fl : 


where f{m,t) is the particle number density per unit 
mass for mass m, and the collisional kernel K{m,m') 
contains all of the relevant information about the inter¬ 
acting masses m and m', such as their mutual relative 
velocities, collisional cross-sections, and sticking efficien¬ 
cies, and can inclu de other properties suc h as particle 
porosities (e.g., see lEstrada fc Cu"^ l2008( l . The diffi¬ 
culty with the application of Eq. (l22)l that has made 
global models of nebula evolution intractable is that the 
calculation is computationally expensive. Coagulative 
grain growth can span many decades of particle mass 
which may require ^ 10 ^ — 10 ^ bins or more depend¬ 
ing on desired accuracy, and at least the first integral of 
Eq. (l22ll must be solved for each mass bin, thus solving 
Eq. (l2^ at every spatial location and at every time step 
can become intractable. It has been shown that taking 
shortcuts with the number of mass bins risks producing 
a rtificial growth. 

lEstrada fc Cuzzil ()2008l 120091 1 developed a scheme for 
modeling coagulative growth that overcomes these diffi¬ 
culties and that lends itself quite nicely to problems in 
which one is interested in globally simulating the evolu¬ 
tion of solids over a wide range of sizes where the de¬ 
tailed size distribution at small sizes both has a gener¬ 
ally predictable form, and is of secondary interest. This 
approach uses a finite number of moments Mi of the 
particle mass distribution, defined by 

poc 

Ml = / m^f{m,t)dm, (23) 

Jo 


dmi^ 

dt 


= {3-q){2-q)pT2 
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^ — m. . ^ 1 
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3-q 


(25) 


For 9 < 2, mn is represented by the upper end of what we 
will refer to as the “dust” mass distribution. For 9 = 2, 
there is equal mass per decade, whereas for 9 > 2 , most of 
the mass is in the smaller particles and the mass mu will 
depend on the smallest size in the distribution. In Eq. 
(Ell), Ta is an integral over the collisional kernel derive d 
from Eqns. (I231I24I) for 1 = 2 f see lEstrada fc Cuzzill2008ll : 


/•mnd) /•mnit) 

T 2 = / K{m,m')m^~'^m'^~‘^ dmdm'. 

(26) 

The d efinition of the kernel K{m,m') is discussed in Sec. 
12.4.21 We solve Eqns. (l2^ and (l26l) using a 4th order 
Runge-Kutta method. 

Because we assume a fixed power law representation of 
the dust particle mass distribution, its mass histogram is 
defined by the moments at any time, the index 9 , and for 
any resolution, by its lower and upper bounds. We define 
a particle radius distribution logarithmically spaced from 
lower bound to upper bound tl: 

fc-i 

'^k — Z'^min) ( 27 ) 


to track general properties of the particle population over 
time under the assumption that the general form of the 
particle mass distribution up to some fragmentation mass 
TO, is known. Note that, by this definition of /(to), the 
mass volume density in a bin of width dm is mf(rn) dm. 
In the moments method, the coagulation equation is re¬ 
duced to a set of ordinary differential equations 


dMi 

dt 





-(to -I- TO ) 



X 


K(m, m') f {m,t) f (m',t) dm dm' 


(24) 


which leads to a closed set of equations as long as the to¬ 
tal number of powers of to in the expression on the RHS 
is < I, which is true even for the realistic collisional kernel 
we use. In this work, we assume that the dust mass dis¬ 
tribution is a power law /(to) oc m~‘^ with (potentially 
variable) exponent q. The assumption of a power law 
distribution for the dust populat ion is motivated by a 
number of detailed models (e-g.. IWeidenschillind 1997, 
20001: iDullcmond fc DominiB 120051 : iBrauer et al.l I2OO8 : 
Birnstiel et al.l 120101 7 120111 12012^ that show distribu¬ 
tions with nearly constant mass per decade up to some 
upper limit which grows with time until a frustration 
limit is reached Isee lEstrada fc Cuzzill2008ll . 


where rin defines the number of point s per decade of ra¬ 
dius. iDrazkowska et aH (120131 120141) have emphasized 
the importance of maintaining good mass resolution in 
brute-force solutions of the coagulation equation, finding 
that ^ 40 bins per mass decade (or ^ 120 per radius 
decade) are generally required. However, such a large 
number of bins is not required for the moments method, 
given our assumption of a power law distribution in the 
dust component. We typically use 20 — 100 bins per ra¬ 
dius decade to calculate the relative velocities, which is 
more than sufficient. The particle masses can then be 
determined using the average dust particle material den¬ 
sity Pp, which we determine from the fractional masses 
of solid species: 


Pp = 




(28) 


where pi is the material density of species i (see Table 2). 
In this paper, we always assume that the initial distribu¬ 
tion has r„iin = 0.1 pm, and rL(t = 0) = 1 pm at all R. 
Though the density of dust particles evolves as their com¬ 
position changes, and thus TOmin is not a true constant, 
the minimum radius r„iin is a constant in our evolutions. 
Much like modeling a variable 9 , we could in the fu¬ 
ture model a different or variable rmin by employing the 
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appro priate number of moments (see lEstrada fc Cuz^ 
l2008f) . The logarithmic binning in (and thus nik) is 
used as the basis for calculations of relative velocities 
(Sec. I2.4.3|l . for characterizing the reservoir of “dust” 
material that migrators may sweep up (Sec. 12.4.41) . and 
to compute the opacity (Sec. 12.3.21) . 

We also note that, it is fairly straightforward to imple¬ 
ment particle_goros^in_the growth parts of our code 
(see lEstrada &: Cuzzil 1200811. Particle porosity may be 
very i mportant (see, e.g. JOrmel et al.ll200^ IZsom et all 
120101 : lOkuzumi et al.l I2012ll in allowing for particles to 
effectively grow to larger sizes while maintaining low St. 
Thus, their radial drift times would be longer perhaps 
allowing them to circumvent the radial drift barrier (see 
Sec. 12.4.51 and 14.11) . Furthermore, the particle poros¬ 
ity can itself evolve with size which would also affect the 
opacity in addition to stopping times. Our code can han¬ 
dle this modification, but we leave this further layer of 
complexity for a later paper. 

2 . 4 . 2 . Particle Sticking, Bouncing and Fragmentation 

The collisional kernel K(rn,m') can be factored into 
three components: 


K{m,m') = a{m,m')AVpp{m,m')S{m,m'), (29) 

where a{m,m') = is the colli¬ 

sional cross section between masses m and m', Kq = 
7r(3/47rpp)^/^, and AEpp(m,TO') is their relative velocity. 
The kernel properties depend on the size distribution of 
particles, the individual particle densities, the total mass 
fraction of solids, and ambient nebula conditions. 

We define sticking coefficients S{m,m'), which depend 
on the particle masses and relative velocities, to capture 
both the “bouncing barrier” and the “fragmentation bar¬ 
rier” as follows. The first barrier encountered by growing 
particles is the bouncing barrier, when S\,{m,m') —>■ 0 
according to: 


Sh{m, m!) = 1 — 


m A’^Vpp(m,m') 
+ m' 


> 0 . 


(30) 


We adopt a b ouncing presc r iption using the second row 
of Fig. 11 of IGiittler et al.l (|2010ll in which the thresh¬ 
old velocity for bouncing collisions between similar-sized 
compact silicate aggregates can be approximated as 14 = 
(Cq/to')^^^, where the constant Cq = 10“^ g cm^ s“^. 
Particle pairs that have = 0 are considered to have 
sufficient energy to avoid sticking, but insufficient energy 
to fragment. 

In a similar way we adopt a condition for the fragmen¬ 
tation of a target with mass m' by a projectile of mass 
m\ 


St,{m, m') = 1 — 


m A^^pp( m, m') 
m m' Q, 


> 0 . 


(31) 


Equations dsni) and (EU account fo r bouncing and 
fragmentatio n criteria implicitly (e.g., iWindmark et al.l 
I20i2al [201311 by smoothly decreasing ;5b,*( m, m') from 1 
at zero collision velocity to 0 for a particle pair colliding 
with a mass-dependent collision velocity. 


The particle fragmentation strength is captured through 
the parameter Q*, which has the dimensions of velocity 
square d as i n Eo. (l30l). We fo llow iStewart fc Leinhardd 
(|2009f) and IBeitz et ^ (1201111 for weak silicate parti¬ 
cles of comparable mass, colliding at low relative veloc¬ 
ity. We include a compositional variation in Q,, moti¬ 
vated by recent results that suggest that icy particles are 
“stickier” or stronger, and might grow larger , faster, and 
with higher porosities than silica te particles (iWada et al.l 
120091 [2^^ lOkuzumi et al.ll2012L see also Sec. [T|). We de¬ 
termine the local value of Q* by a mass weighted average 
over the species making up the composition of a grain 




(32) 


Because we consider icy particles to be stickier, we also 
scale the bouncing threshold velocity by a factor of 10 
(or 100 in energy) as we do for the fragmentation. Al- 
though the kernel can a lso include particle porosities (see 
lEstrada fc Cuzzil 1200811 . we do not include them explic¬ 
itly for this paper. In practice this leads to Q* — 4 x 10® 
outside the ice line, and Q* = 10"^ inside. 

The fragmentation barrier mass to* is determined by 
the condition that the energy per unit mass in a collision 
between a target particle of mass to' = to* and some 
other particle to < to* exceeds the strength of the parti¬ 
cle Q*. In turbulent conditions, to* is the mass of a par¬ 
ticle that is destroyed upon colliding with a comparable 
mass particle. However, the radius of the particle that 
satisfies the condition in Eq. may be smaller than 
TO* under low turbulence conditions where headwind- 
drag-driven velocities dominate (See next section). 

Once the fragmentation barrier is reached for a tar¬ 
get particle, our code provides a reservoir of material 
(with an associated creation rate) from which growth 
may proceed by sweep up of smaller grains, but which 
remains subject to destruction by particles of compa¬ 
rable or smaller sizes. These particles represent the 
lower end of the submigrator population. Particles that 
are fragmented release their mass into the background 
“dust” population with its current local size distribution. 
We treat the fragmentation of particles statistically (see 
Sec. I2.4.'B1) as well as mass transfer between them (Sec. 
12.4.41) . Thus we believe our model captures the essen¬ 
tial physics of recent experimental outcomes and models 


tion (IGiittler et al.1 

2?i09l [20inl: IZsom et al.l 

2nTnl IWll: 

IWeidling et al.ll2012 

; IWindmark et al.ll2012a 

). 


2 . 4 . 3 . Relative Partiele Veloeities 

We include a variety of sources for the particle rela¬ 
tive velocities: Brownian motion, pressure gradient, ver¬ 
tical settling and turbulence. We briefly summarize them 
here, while giving a more detailed description of how we 
calculate them in Appendix A.4. 

The thermal motion of particles. Brownian motion, is 
dependent upon the masses of the particles to and to', 
and the ambient nebula temperature 


AVBro(in, to') 


SfeeT (to + to') 


mm' 


(33) 


and is only effective for the smaller particles near the 
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lower bound of our mass distribution. 

The pressure-induced, or systematic dust velocities re¬ 
sult from the gas in the nebula orbiting at slightly less 
than the local Kepler velocity. In a rotating frame, a 
parcel of gas experiences an outward directed pressure 
gradient force that counters the inward force of solar 
gravity. The dust particles do not feel this radial pres¬ 
sure gradient force directly, but experience an azimuthal 
drag force from the more slowly rotating gas, leading 
to size-dependent radial and azimuthal velocities. In 
cases where the local solids fraction is high and affects 
the gas velocity, particle and gas relative velocities must 
be determined iteratively. To ensure correct results in 
all regimes, we routinely solve for these components of 
the velocity from a s et o f equations generalized from 
iNakagawa et al.l (jI98fifl bv lEstrada fc Cu^ (120081) for a 
parti cle size distribution Ithough. also see iTanaka et al.l 
[ 2 ?) 0 ^ : 

^ = -Akp{Uk-u) + 2nVk ; (34) 

at 

dVk I 

= -Akp{Vk - v) - -ftUk, (35) 

^ = + (36) 


Uk 


2Stkr]VK 

1 + Stfc 


(40) 


For the turbulence-induced velocities, we use the closed 
form prescripti o ns fo r a particle size distribution of 
lOrmel fc Cuzzil (I2007f) . The relative velocity with re¬ 
spect to the gas is th en (|Cuzzi fc HoganI 

120031 : lOrmel fc Cuzzill200';^ . where the turbulent gas ve- 

locity is given by Vt = (see, e.g.. lCuzzi et al.l[200I[l . 
and Vt is the average inertial space particle velocity due 
to turbulence. The particl e-to-particle turbulen t relative 
velocities AV 12 (Eq. ('161 lOrmiel fc Cu^l2007|l are less 
straightforward because of the different coupling that ex¬ 
ists between particles and eddies of different sizes. Re¬ 
cent numerical simulations have obtained results differing 
from these para metrizations by a factor of order unity de¬ 
pending on St (|Hubbardll 20 I^ [Pan fc Padoanll20I3ll . It 
remains unclear how much of this difference arises from 
the relatively low Reynolds number of the numerical sim¬ 
ulations, relative to the inertial range turbulent kinetic 
energy spectrum of the actual nebula which is assumed 
by Ormel and Cuzzi (2007). For more discussion, see 
iCuzzi &: HoganI (1200311 . 

From the various contributions, we calculate the 
particle-to-gas relative velocities for a particle of mass 
m which we use to calculate the particle stopping times; 


^ = -^Aip^(v-V0-^nu, (37) 

where (Uk, I 4 ) are the radial and azimuthal velocity com¬ 
ponents for particles of mass mk, (u, v) are those for the 
gas, Ak = (ptl)~^, and 


AEpg = [(U - uf + {V- vf + 1/2 ^ 

where W = —Vt^zts is the particle vertical settling veloc¬ 
ity, in which the vertical coordinate z is generally taken 
to be the subdisk scale height Hd (see Sec. I2.2.3|) . The 
particle-to-particle relative velocities are then given by 


Pk = 


V TT 2/1^ 


(38) 


is the mass volume density of mass bin k where is the 
scale height of ruk (see Sec. 12.2.31) . 

Equations (I34II37|) represent a set of 2n -I- 2 equations 
in 2 n -|- 2 unknowns where n is the number of particle 
bins in the distribution given that there are rip parti¬ 
cle bins per decade radius (see Sec. 12.4.II and Appendix 
B). We solve this system of equations using a matrix 
method as defined in Appendix A.4. The pressure gra¬ 
dient in Eq. (I36L where p is the gas pressure, can be 
expressed in terms of the more familiar p parameter as 
— (1/p)dp/dR = 2II?7Ek, where 


p(R,t) 


1 

2 


/ c y q 

VEcj Eri/2 ^ 


^3/2j.j.l/2 


(39) 


and Vk is the local K epler velocity (|Nakagawa et al.l 
II986t[Ouzzi et al.llI993ll . Pressure gradients can be quite 
steep near the outer edge of the disk, and can lead to 
rapid inward migration of even very small particles. In 
most cases, when the local particle density is small com¬ 
pared to the gas density, the particle drift velocity Uk can 
be well approximated by (Weidenschilling 1977, Cuzzi 
and Weidenschilling 2006): 


AVpp(m,m') 


(AEBro)" + (APpre)" + (AE12)" 


(42) 

where (AVp^ef = (U - U'f + (V - V'f + (W - W'f. 
If we are strictly in the Epstein flow regime, then the 
stopping times do not depend on AEpg, so AEpp can be 
calculated in a straightforward manner from the velocity 
components. However, if there are particles in the Stokes 
flow regime, then the stopping times depend on AVpg and 
iterations are needed to converge to a proper solution (see 
Appendix A.4). 

We note that in Eq. m, we have summed the mean 
relative velocity between the gas and particle in quadra¬ 
ture with that of the fluctuating turbulent relative ve¬ 
locity which is not strictly correct. A similar argument 
can be made for the inclusion of AVpre in Eq. (l42)l which 
we take t o be the mean collis ion velocity. It has been 
argued bv iGaraud et al.l (j20I3ll that this construct is not 
accurate when the dominant velocities are systematic. 
However, in the models we present here, the turbulent 
velocities dominate the systematic ones so that we do 
not expect that any differences will be significant. We 
leave a more proper treatment for future work. 


2.4.4. Migrators and Growth Beyond the Fragmentation 
Barrier by Mass Transfer 
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Once the fragmentation barrier m* has been reached, 
we employ a more sophisticated, semi-analytical and sta¬ 
tistical algorithm for the further growth of particles. In 
this regime, we track mass and radial drift explicitly for a 
population ranging from “submigrators” which are sub¬ 
ject to mutual collisional destruction from particles of 
similar size or smaller, to “migrators”, which have grown 
large enough to have a much lower probability of de¬ 
struction and may grow by ma ss transfer between them 
and particles of s maller size (jWindmark et all l2012al 
iGaraud et al.ll2013[) . Disrupted submigrators and migra¬ 
tors are distributed back into the dust population with 
the same q. The migrator distribution is defined on a 
continuous, time-dependent mass grid where the width 
of each mass bin can vary due to different rates of growth 
for different sized particles. Thus the evolving mass his¬ 
togram for m > TO, does not follow a power law in gen¬ 
eral. 

Beyond the fragmentation barrier, particles can still 
grow through incremental means, by the sweepup of 
smaller material. The simplest expression for incremen¬ 
tal growth assumes perfect sticking between a particle 
and (smaller) feedstock particles so that, schematically, 

^ = aAVdPd, (43) 

where tr is the collision cross section and AVd is 
some characteristic relative velocity between to and 
the feedstock pop ulation (see, e.g., iCuzzi et'all 119931 : 
iBrauer et al.ll200§l . We generalize this to a size distri¬ 
bution, in which sticking is not assumed to be perfect, to 
define the growth rate of a migrator mk > : 


drrik 

dt 


iL 

Y,<yk,iSk,iAV^!pf 


k 

+ '^k,iSk,iAVk^^pT, 

Z=n* + 1 


where Sk,i are the sticking coefficients, the index n* refers 
to the fragmentation barrier mass bin, and AV^^ is the 
relative velocity between ruk and to;. The first term on 
the RHS is due to the dust population, and the second 
term the migrator population. The volume density of 
migrators p™ is defined in terms of the total mass of 
solids in a migrator mass bin, and the total volume 
of the particle sublayer, of vertical thickness 2hd near 
the midplane, from where migrators can accrete other 
material 


Pk — 


2^/hd ’ 


(45) 


where is the surface area of the radial bin. Equation 
allows for a migrator to^ to accrete other particles of 
"fni < mk, but the sticking will be zero for a large range 
of size pairs when the collision specific energy exceeds Q, 
(or Co). 

On the other hand, there are outcomes of high-velocity 
collisions which lead to growth of the target particle, 
which have been observ e d and studied exp erimentally 
(e.g., IWurm et al.l 120051 : iKothe et al.l 1201011 . Specifi¬ 
cally, the impact velocity may be insufficient to frag¬ 
ment the larger target particle, but sufficient to fragment 
the smaller particle and allow deposition of mass on the 


larger particle - if the efficiency of accr etion e„r. exceeds 
that of erosion Cer. We use the model of lWindmark et aP 
dMHB), in which the threshold for fragmentation is 
both mass and velocity dependent, to account for this 
process: 


M(m, Wm) = (46) 

where Ecm is the relative velocity of the target and 
impactor in the center of mass frame, the coefficient 
Cw = 3.27 g- 0.068 masses and velocities are in cgs 
units. In order to account for the stickiness of icy parti¬ 
cles (see Sec. 12.4.21) . we scale Cw by a factor (Q*/10^)°-^^. 
The fragmentation threshold is reached when /i = 1. We 
apply Eq. (l46l) to both to/ and mk where the center of 
mass velocities are given by 

; Vr ^ ■ (47) 

1-l-TOfc/TO/ l-|-TO//TOfe 

Fragmentation with mas s transfer only occurs w hen 
Pfc > 1 and pi < 1 (see iWindmark et al.l 1201^ . If 
this condition is sati sfied, we then calc u late the efficiency 
of acc retion from (|Beitz et al.l 120111 : IWindmark et al.l 
I2012all 

/ xO.16 

Cac = -6.8 X 10-3 + 2.8 X 10-4 ( ^ ) (48) 

and the efficiency of erosion (IWindmark et al1l2012af) 

/ ^0.15 

Cer = 9.3 X 10-3 — AEpp - 0.4, (49) 

\moJ 

where toq = 3.5 x 10-^^ g is a monomer mass. The 
erosion efficiency is an interpolatio n between the exper¬ 
imental results of several workers ([Paraskov et al l 120071 : 
iTeiser fc WurmI 120091: iSchraoler fc Wurmll20lH) . A suc¬ 
cessful transfer of mass to the larger particle mk occurs 
then if Amk/mi = Cac — Cer > 0, which is assigned to the 
value of Sk,i in Eq. (|4^ . Inspection of Eqns. (l48l) and 
Eq. demonstrates that for a given impactor mass to/ , 
higher relative velocities are required to initiate erosion 
versus accretion. As an example, for a mk = 1 g target 
particle being impacted by a projectile with to./ = 0.05 
g for a relative velocity of 5 m s-^, the target particle 
accretes roughly 6% of the impactor with no net erosion. 
On the other hand, for the same pair impacting at 15 
m s-^, the target particle accretes over 20% of the im¬ 
pactor, but also experiences ^ 7% erosion. In our treat¬ 
ment, the remaining projectile mass, which by definition 
is fragmented, as well as the eroded mass are assumed to 
be returned to the dust population. 

Finally, once migrators are large enough that A^Vpp < 
2(5*, they can accrete migrators as large as themselves in 
pairwise mergers if they collide. A check is made to en¬ 
sure that no more mass is “accreted” in a timestep than 
actually exists locally (for more detail, see Appendix 
A.5). 


2.4.5. Radial Drift and Evaporation Fronts (EEs) 
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Due to variable coupling with the nebula gas, parti¬ 
cles of different size will drift radially at different rates 
with respect to the gas. The radial drift velocity for 
a particle of mass ruk. has two contributi ons (e.g., see 
iTakeuchi fc Linl[2002t iBirnstiel et al.l[2M^ : 

= Y^+At/fe. (50) 

The first term is directly imposed by the radial motion 
of t he g as that moves with advective velocity Vg (sec¬ 
tion IQ) . The second term {AUk = Uk — u < 0) is the 
radial drift velocity of the particle with respect to the 
gas (section [2A3]) . This drift increases with particle size 
until ts ~ or St = tsD ^ 1, but then decreases 
for larger sizes (see Eq. [40]). In the terrestrial planet 
region, meter size particles drift most rapidly, but fur¬ 
ther out in the disk where gas densities are low and the 
pressure gradient can be very strong , much smaller par¬ 
ticles drift inward the mo st rapidly (iBrauer et al.l [20081 : 
iHughes &: Armitag3l2010ll . 

In our code we track the radial drift across radial bin 
boundaries for individual migrator mass bins. We do 
this by calculating the inward drift time of the migrator 
across (logarithmically spaced) radial bins of width AR: 

= ( 51 ) 

where the radial drift velocity for larger particles is gen¬ 
erally a negative quantity (consistent with our sign con¬ 
vention, see Sec. l2.Il or [^7^^ . This allows us to calculate 
for every migrator mass bin k how much mass is drifting 
out of the local radial bin in a time At{R): 

= (52) 

where the factor of 2 assumes that the drift time from 
the bin’s midpoint is representative. Equation (1521) is 
the maximum amount of mass that can drift out of a 
radial bin in time At, and does not take into account the 
fragmentation of some particles (see Sec. 12.4.611 . 

It is possible that some migrators with m ^ m* have 
positive outward drift if, for example, at is sufficiently 
large (which leads to more vigorous collisions and smaller 
TO*) that the gas advection term (first term on the RHS 
of Eq. [50] 1 can overwhelm the radial component of 
the headwind-driven drift velocity. Under such circum¬ 
stances, the redistribution of these migrators is done us¬ 
ing the radial diffusion-advection equation (Sec. 12 .2.21) . 

Self-consistently modeling a globally evolving nebula 
requires that we consider the presence of Evaporation 
Fronts (EEs) which are regions in the disk where phase 
changes between solids and vapor can occur. The evapo¬ 
ration of radially drifting material, and subsequent recon¬ 
densation of outwardly diffusing vapor, can have three 
main effects. First, it can increase the abundance of va¬ 
por inside an EF, with implications for chemistry and 
mineralogy. Second, it can increase the fractional abun¬ 
dance of solid material available just out side the EF, per¬ 
haps by a factor of ^ 10 or more (e.g.. iCuzzi fc Zaimi^ 
l2004t iCiesla fc Cuzzil[20^ lGaraudll2007D . with imolica- 
tions for accretion. Third, it can significantly change 
the composition of solid material outside the EF from 


“cosmic abundance” values. We will illustrate all these 
effects in this paper; a more detailed study will await 
future publications. 

Rather than treating an EF as a sharp boundary, we 
allow for evaporation (or condensation) to occur over 
a small range of radii covering a midplane temperature 
range 2 ATef O.I — I K relative to the nominal species 
evaporation temperature Ti (see Table 2). The fractional 
solid abundance a''’^ of a species i (with density pi) at 
some radial location R and at some local temperature T, 
is transformed into vapor linearly as the calculated mid¬ 
plane temperature changes over the temperature range 
Ti — ATef at < Ti ATef (Sec. 12. 3. II) . This grad¬ 
ual radial transition in solids abundance, which can span 
a significant radial range, mimics the anticipated effect 
(as described below) of “buffered” temperature changes 
across an EF as material is evaporated or condensed, 
rather than allowing unrealistic (and numerically prob¬ 
lematic) abrupt radial changes in opacity and tempera¬ 
ture just inside an EF. This simple numerical treatment 
captures the essence of the actual condensation process 
in which material first evaporates at the midplane, and 
then at increasingly higher altitudes with de creasing dis¬ 
tance from the s tar as the disk gets warmer (|Davisll2005l : 
iMin et aI1l20IIh . The effect is seen in the constant mid¬ 
plane temperature regions in our simulations (e.g., see 
Fig. (Tj). 

In our code, dust grains are effectively treated as ag¬ 
gregates of chemically distinct monomers (although see 
Sec. 12.4. II) whose fraction of species i can be quickly re¬ 
moved or emplaced. Larger particles that are followed 
explicitly (Sec. 12.4.41) are assumed to lose that fraction 
of their material that is of evaporating or condensing 
species i as they migrate through the EF for species i, 
and their masses and mean densities are adjusted accord¬ 
ingly. Their evaporated material is then added to the 
local vapor inventory. Better models of the largest parti¬ 
cles whose interiors are somewhat insulated from ambient 
nebular conditions, and could potentially transport very 
vol atile species from v ery cold to warmer regions (e.g., 
see lEstrada et al.l[2009H , will require physical evaporation 
rates and internal structure models. We do not treat the 
kinetics of evaporation and c ondensation here (se e dis¬ 
cussions in [CuzireFaI][200l; [Cieil£(^^uii^|20^. We 
plan to include these effects in future work. 


2 . 4 . 6 . Probability of Destruction and “Lucky Particles” 

In our treatment of growth, the particle mass such that 
the sticking coefficient first approaches zero for equal 
mass particles represents the “bouncing” barrier. That 
is, such a collision is energetic enough to prevent any 
sort of sticking, but not energetic enough to fragmen t 
the particle (e.g., iGiittler et al. I2010l:[^om et al.ll2010ll . 
Our nominal fragmentation mass to* (section 12.4.21) is 
also a convenient reference value, at which particles of 
equal mass fragment each other under local nebula con¬ 
ditions. In practice, we employ a statistical scheme for 
the probability of destruction of a migrator. The scheme 
utilizes a Gaussian PDF of relative velocities for a given 
impactor mass, with rms velocity equal to the mean 
([Garballido et al.l I20I0I: iHubbardI I20I2I : iPan fc PadoanI 
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nnn). 

The collision rate of target migrators of mass m! with 
projectile particles of masses m < m' is given by {cf. Eq. 

m) 

, m) = 7r(/ + r)^b{m)AVpp{m', m), (53) 

where b{m) = pd,ml'm is the number density of particles 
of mass m, and the subscripts refer to either the volume 
density of dust or migrators. The probability that a mi¬ 
grator of mass m' will suffer fragmentation as it grows 
during time At can then be obtained from 

nm' 

g^{m')= / dt' 7r{r'+ r)'^AVpp{m',m) 

Jt Jmo ’ 

xC(m', m)f{m, t') dm 

where (^{m',m) is an integral over a Gaussian dis¬ 
tribution of relative velocities covering a range equal 
to or greater than the critical impact velocity 14 = 
V2Q ,(to + m')lm: 



-1 poo (V'-AVpp)^ 

C(m',m)= - - e dV 

V27rAEpp(m',m) Jvcim) 

14-AEpp V ■ 

V2AVpp j 

(55) 


1 — erf 


We discretize Eq. dSl for use in our calculations 
as described in Appendix A.6. Although other work¬ 
ers have utilized a Gaussian scheme as we do here, 
others have arg ued that the distribution should be 
Maxw elli an (e.g., iGalvagni et al .1120111 : iWindmark et al.l 
l2012bll . iGaraud et al.l ( 201311 have developed a more 
detailed model in which she shows that the mean and 
rms square velocities are not the same. Furthermore, 
our approach to calcula ting C, is more akin to t hat o f 
IWindmark et al.l ()2012b[ l. whereas IGaraud et al.l ()2ni3[ l 
argues that the relative velocity (in our notation V') 
should be included in the integral over the PDF in order 
to recognize that collisions of particles with larger col¬ 
lision velocities have a higher collisional frequency than 
those with lower collision velocities. We intend to exam¬ 
ine the more detailed model of Garaud et al. in future 
work. 

Using our formalism, we calculate the fraction of mi¬ 
grators within every migrator mass bin that is destroyed 
during every growth step. The amount of mass re¬ 
turned to the dust distribution for any mass m' is then 
Afdest = ^{m’)M^{m') where is the total mass 

in migrators with mass m'. This mass is subtracted from 
the total mass contained within the m' mass bin prior to 
determining how much mass drifts out of a radial bin 
(and thus a factor 1 — is included in Eq. [5^1. 

We find quite generally that growth stalls at masses 
slightly-to-moderately larger than the fragmentation bar¬ 
rier mass m* and radius r*. We will represent this par¬ 
ticle radius where growth stalls to be tm which also is 
the particle containing the bulk of the total mass of mi¬ 
grators in a radial bin. This result of stalled growth is 


Figure 1. Initial profiles of the Rosseland mean opacity (red 
curve), and temperatures at the midplane (solid black curve) and 
photosphere (dashed curve) for our fiducial model with at = 
4 X 10“^. Additional midplane temperature curves are shown (dot¬ 
ted curves) for different values of at- The nearly flat regions in the 
temperature profiles correspond to the location of EFs as labeled. 


fairly robust unless the nebula turbulence is very small, 
in which case incremental growth can proceed to large 
sizes because the particle layer becomes very dense near 
the midplane, and coll ision velocities ar e low. 

As was seen by IWindmark et al.l (l2012al lbl) and 
IGaraud et al.l (120131 ). we also find that a small fraction 
of migrators are able to grow to larger sizes because they 
have grown large enough to have a low probability of 
destruction. We refer to these migrators as “lucky parti¬ 
cles” , and characterize the largest particle in each radial 
bin by radius rp. However, previous studies (except for 
Drazkowska et al 2013, who imposed a “pressure bump” 
to prevent loss by drift) have not allowed for the signifi¬ 
cant radial drift of such particles, which greatly reduces 
their abundance. While rare, such particles can trans¬ 
port material over large radial distances. Treating the 
growth of these “lucky” migrating particles simultane¬ 
ously with their radial drift and nebula evolution is a 
significant advantage of our code over previous global 
models. 


3. RESULTS 
3.1. Baseline model 

We choose as our fiducial model a nebula with a stellar 
mass of M* = 1 Mq, an initial disk mass of 0.2 Mq and 
a scaling factor Rq = 10 AU (see Sec. lATD . The radial 
grid in our evolutions spans a range from 0.5 — 1000 AU 
using 96 logarithmically spaced bins, chosen to make op¬ 
timal use of the Haswell CPUs (12-core, 24 processors) at 
the NAS Pleiades cluster used in our simulations. In our 
code, radial bin width is only a concern for radial drift 
of migrator particles, because other properties that com¬ 
municate across the grid such as gas and dust evolution 
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are fully implicifl All simulations in this paper span up 
to 2 X 10^ years, and use variable stellar luminosity L 
and particle strength Q* (except in sections 13.51 and 13.61 
where these parameters are kept constant). 

We employ a constant turbulent viscosity parameter 
at = 4 X 10“^ as our baseline, though we explore a model 
that uses a constant at = 4 x 10“^. As noted in Section 
[U recent studies of purely hydrodynamical instabilities 
suggest that at > 10“"^ (Nelson et al 2013, Marcus et 
al 2013, 2015; Stoll and Kley 2014). We do not model 
disks with lower value s of at. Such models have been 
explo red in detail le.g.. iWeidenschillindllOOTl 12000112004 
Mil : vanishingly small at leads to rapid in situ growth 
by sticking and/or collective instabilities rather than the 
delayed growth, and extensive redistribution of solids, 
which is the primary focus for this paper and arguably 
in better agreement with observations (Sec. [I}. 

In the following sections, we consider models with frag¬ 
mentation (F) only, then bouncing and fragmentation 
(BF) and finally add mass transfer (MTBF) as described 
in Sections 12.4.21 and 12.4.41 In forthcoming papers, we 
will explore radially and vertically variable at profiles 
le.g.. IBai fc Stonjl20lH iKalvaan fc Deschl[201^ . which 
introduce another layer of complexity. 

In this work, we choose to cut off the initial distri¬ 
bution of solids at 100 AU (where the initial gas surface 
density E quickly drops off to <C 1 g cm~^). Other work¬ 
ers have chosen sim ilar cutoffs (e.g., iBrauer et al.l 120081 : 
iBirnstiel et al.l [201011 : while arbitrary, this initial condi¬ 
tion simply provides a convenient reference point from 
which to follow the solids evolution. The motivation is 
that the strong pressure gradients at the outer edge of 
the nebula will lead to rapid inward migration of even 
the smallest grains unless one assumes sufficiently high 
values of at which leads to outward transport of solids as 
the gas disk spreads in spite of the presence of a strong 
pressure gradient (e.g., see Sec. 1X71) . 

In Figure [T] we show the initial midplane temperature 
and opacity profile for our fiducial model (solid curves). 
The corresponding photosphere temperature is shown by 
the dashed curve. The initial disk temperature is suffi¬ 
ciently hot at the midplane for our high at model (Sec. 
13.71) that all of the EFs for the species used in this paper 
(see Table 2) are present. Although midplane tempera¬ 
tures in EF regions appear flat, there is a small gradi¬ 
ent in temperature (Sec. 12.4.51) over which evaporation 
of volatiles is assumed to take place. It is over these 
transition radial ranges that the opacity decreases lin¬ 
early. For comparison, we plot the temperature profiles 
for different values of at (dotted lines) to the show the 
sensitivity to the choice of this parameter. Even for a 
value of at = I0“®, the disk is relatively hot in the inner 
regions because the initial ste llar luminosity is ^ 12 Lq 
([D Antona fc Mazzitellil[l994[) . In this work we adopt a 
constant value of 10“'^ cm^ g“^ for the gaseous opacity. 
In the hot, innermost regions this may be an underesti¬ 
mate beca use gaseous opacity i s a function of tempera¬ 
ture fe.g.. iFerguson et '^11200511 . However, for our mod¬ 
els, dust is usually present at all locations in the disk 

^ A minimum requirement we impose is that the timestep used 
in our simulations is much smaller than the drift time of the fastest 
migrating particle across a radial bin. For all of our simulations 
presented here, this is easily satisfied. 


t = 2x10^ yeors 





Figure 2. This figure shows several of the actual histograms of 
particle size distributions, as mass volume density per bin of width 
dm, for three different cases and several distances from the star. 
The lower boundary of these plots aligns with the upper boundary 
of the “dust” power law distribution at r*, that extends downward 
to O.l/rm radius. The metho d of determining these distributions 
is described in section 13.1.11 As is specifically indicated in the 
upper panel, the peak in this distribution represents the radius 
rM which characterizes most of the net mass. The largest single 
particle in each distribution, at the far right end, provides our value 
of • It is clear that the mass contained in these “lucky” particles 
is extremely small. In the lower panel one sees secondary peaks 
in ( g cm~^), which a re like ly du e to mass transfer a s 

suggested by IWindmark et al.l Il2012al^ and IGaraud et aT\ 1I2013IV 
Though we might regard these as “breakthrough” particles, they 
are actually not that much larger than tm , and still contain much 
less mass than particles of size rjvi- This figure should be kept in 
mind for interpreting the remaining figures where only rM and rL 
are shown. 

SO that the Rosseland mean opacity is dominated by this 
component. In future work, we will incorporate a gaseous 
opacity table into our models. 

3.1.1. Typical particle size distributions 

Before discussing the different cases in detail, we show 
some typical particle size distributions in Figure 
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Three of our cases (discussed below) are shown in the 
three panels. The curves in each panel are actual calcu¬ 
lated distributions of particle mass volume density as a 
function of particle radius for particles with r > r*, cal¬ 
culated with 100 bins per decade radius; recall that par¬ 
ticles with r < r* are modeled by a “dust” distribution 
which is a powerlaw extending downward and to the left 
of the values shown. Each family of curves is for several 
different subsolar distances between 0.5 — 1 AU. Values 
of the fragmentation radius r,, the mass-containing ra¬ 
dius tm (shown for the red curve only), and the largest 
“lucky” particle radius tl are indicated in the top panel. 
The downward arrow indicates that the largest particle is 
actually below the range plotted which covers seventeen 
orders of magnitude in mass volume density. In future 
plots we will only show these three characteristic values, 
but the distributions of Fig. [2] underlie them all. In the 
lower panel, secondary peaks are seen which are most 
likely due to mass transfer as will be discussed later. 

3.2. Fiducial model with fragmentation (“F”) 

In our model, the fragmentation size is defined by that 
particle mass m' = m* that gives S'* = 0 in Eq. (EH). 
The value of to* will depend on ambient nebular condi¬ 
tions, as well as the particle strength Q*, which is compo¬ 
sition dependent (Eq. E2] ) for the bulk of the simulations 
we present here. 

Figures [3] and [4] show the results of a fragmentation- 
only model ”F” in our nominal, moderately turbulent 
disk. The temperature (black curves) and Rosseland 
mean opacity (red curves) are plotted in the top panel of 
Fig. Oat evolution times of 10^, 10® and 2 x 10® years. 
Early on, the disk remains relatively hot in the inner re¬ 
gions as seen by the persistence of the silicate EF at 1 AU 
even after 10® years, but the temperature quickly drops 
to much lower values by 2 x 10® years. This cooling is a 
direct consequence of the lower opacity associated with 
rapid particle growth when only limited by fragmenta¬ 
tion. Evaporation fronts can clearly be seen to evolve 
over time, with the organics and water ice EFs at 425 
K and 160 K, respectively, being the most prominent at 
later times (solid curve). 

Fractional masses of different components are shown 
in the lower panel of Fig. [31 Here, the red curves refer to 
the dust, magenta to the migrator and green to the va¬ 
por fractions, respectively. Peaks appear in the fractional 
masses of the dust and migrator populations outside of 
EFs, because migrating particles passing through them 
release their volatiles into the vapor component, and 
some of this vapor subsequently diffuses outward across 
the EF to recondense onto grain s urfaces. This effect has 
been seen in previous simulations (iCuzzi fc Zahnlell2004 
iCiesla fc Cuzzil 120061 iGaraudI 120071 ). Relatively strong 
peaks in the fractional masses of solids at 10® years can be 
associated with the organics and water ice EFs, whereas 
at 2 X 10® years these peaks have migrated inwards con¬ 
siderably as the disk cools. A less prominent silicate peak 
is also seen at 10® years, that has all but vanished at later 
times, but this is because the silicate EF is no longer on 
our computational grid, and the inner boundary cannot 
relay information about what is occurring inside itself. 
The organics enhancemeni0 is then largest at around 1.5 

® For this introductory work we assume that “organics” are the 


AU, while the water ice peak lies near 5 AU. The solid 
phase mass enhancements map well with the the opacity 
shown in the top panel. That is, the highest opacities 
in the inner disk at all times are associated with the 
peaks in solids. A steady transport of material from the 
outer regions to the inner regions maintains a significant 
population of smaller particles, which keeps the opacity 
high even as particles grow, and maintains a high tem¬ 
perature. This can be seen by noting that at 10^ years, 
there is still a lot of material in the outer disk out to 100 
AU (large opacities and fractional masses there). By 10® 
years, for these specific models, the bulk of the material 
in the outer nebula has effectively been transported into 
the inner nebula, though some material still remains in 
the outer disk (at fractional masses of < 10“® — 10“®). 
We see overall enhancements in the inner regions of a fac¬ 
tor of ^ 2 — 3, though it is mostly present in the vapor 
component (green curves). The inner nebula situation 
changes significantly by 2 x 10® years where much of the 
material at R < 1 AU has drifted inwards leaving much 
lower fractional masses of solid material, and thus much 
lower opacity leading to cooler overall temperatures. 

The relatively rapid transport of solids to the inner re¬ 
gions occurs because a fragmentation-only model allows 
larger particle sizes to be reached relatively quickly in 
the outer nebula, and in relatively large numbers. In the 
middle panel of Fig. |3| we plot particle sizes of interest 
(black curves), particle Stokes numbers for tm and tl 
( blue curves) and the mean particle density of the dust 
(red) and migrators (magenta) at 2 x 10® years. The first 
thing to notice is that the fragmentation barrier size r* 
has been reached almost everywhere in the disk. This is 
indicated by the presence of a heavy dashed curve at all 
distances (actually the situation was reached at around 
5 X 10^ years). Migrator particles also form early on, 
which accelerates inward solids transport. Once the frag¬ 
mentation barrier has been reached, further incremental 
growth is hampered by the gradual destruction of parti¬ 
cles r > r* with time, and the size of particle that carries 
most of the mass in the migrator size distribution (tm) 
does not stray far from the fragmentation size. The char¬ 
acteristic “hump” outside of 3 — 4 AU corresponds to the 
transition between weaker and less sticky silicate parti¬ 
cles inside the water ice EF, with the stronger and stick¬ 
ier icy particles (with larger Q*) outside the EF. Indeed 
for this case, outside the ice EF, most of the mass has 
reached tm ~ meter size. Although incremental growth 
by sweepup of particles too small to fragment large tar¬ 
gets leads to “lucky” particles approaching tl ~ meter 
size even for some distance inside the ice EF, these par¬ 
ticles are very few in number and contain a negligibly 
small portion of the mass (see Fig. jj). This distribution 
of sizes essentially stays the same between 5 x 10^ and 
2 X 10® years, although there is some variation in transi- 

moderately refractory, carbon-rich materials seen in comet Halley 
and p rimitive carbonaceous chondrites, which vo latilize at roughly 
400K IjGradv &; Wrigh^l2003l:[Pollack et alJ1994l . In what is prob¬ 
ably an oversimplification, we further assume that these materials 
behave like normal volatiles, able to reversibly evaporate and re¬ 
condense. In reality, their “evaporation” might be closer to an irre¬ 
versible decomposition into simpler molecules such as CO and/or 
CH 4 . This refinement is deferred to future work, but there will 
be implications for opacity, solids mass enhancement and nebula 
chemistry. 
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F: vor Q.. cXj = 4x10-* 
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Figure 3. Fiducial model with fragmentation only (F). 
upper panel: Temperature (black) and opacity (red) plotted at 
three different times. middle panel: Fragmentation radius r* 
(dashed), particle radius that carries most of the mass tm (solid 
black) and largest particle radius (solid black curve). The Stokes 
numbers for “Tm and are shown in blue. Also plotted are the 
mean particle densities for dust (red) and migrators (magenta), 
lower panel: Fractional masses shown at two different times. The 
F model is characterized by rapid growth and drift, leading to de¬ 
pletion of the outer disk and only modest enhancements in the 
inner nebula. However, the largest particle sizes can be > 1 m. 
See text for details. 

tion boundaries due to the radial motion of the EFs. 

The particle Stokes numbers (blue curves in the middle 
panel of Fig. [SD help to demonstrate why there is sys¬ 
tematic mass transfer from the outer nebula to the inner 
nebula. The solid Stokes number curve refers to vm and 
the dotted Stokes curve refers to the largest particle in 
the size distribution which include the “lucky” particles 
> ^M- They diverge inside the ice line, where only 
rare lucky particles have such large values of St and tl- 
Even while, near the ice EF, particles of about 1 m radius 
can be found, the high gas densities mean that St < 1 ev¬ 
erywhere except at the outermost edge of the solids disk 


Figure 4. Fiducial model with fragmentation only (F). 
upper panel: Gas (black) and dust (red) surface densities at three 
different times, middle panel: Radial velocities of the gas (solid), 
largest particle (dotted), and mean inward (long dashed) and out¬ 
ward (short dashed) components of the dust population at 2 x 10^ 
years, lower panel: Mass accretion rates for the gas (solid black) 
and the dust component (red curve). The vertical scale for the dust 
component is enhanced by a factor of ten for clarity. The dashed 
curve denotes the pressure gradient. For the fiducial choice of at, 
the evolution of the gas is modest, but substantial changes in the 
solids arise due to rapid growth and inward migration. See text for 
details. 

200 AU) where the gas density is very small. Recall 
from Eq. (gOl) that for St < 1, the radial drift velocity 
is proportional to St. Inside the ice EF, StM drops by 
almost two orders of magnitude because of the smaller 
Q*, causing the bulk of particles to drift inward more 
slowly than outside the ice EF. This explains why there 
is a pileup of solid material inside the ice EF. The fact 
that tl gradually decreases from just outside the water 
ice EF to the inner nebula is due to the fact that even the 
lucky particles are not growing fast enough to overcome 
the radial drift barrier. This is discussed in more detail 
in Sec. |4](see Fig. 1181) . The maximum drift rate for par- 
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tides is at St = 1 (Eq. [ID]), but particles with St < 1 are 
still drifting faster as they grow larger, so these particles 
are suffering a gradual but steady loss inwards. Because 
the radial drift barrier cannot be overcome means that 
eventually most all material will be lost at later times. 

Interestingly, the maximum Stokes number Stn is 
roughly constant outside the water ice EF, at around 
a value of ^ 0.1. A simple relationship explains 
this behavior in the fr agmentation-dominated r egime 
(IWeidenschillin d Il98^ iCuzzi &: Weidenschillind 120061 
iBirnstiel et al.l l2009l l201l! h One merely sets the col¬ 
lision speed AEpp equal to the fragmentation speed 
Uf ^ Since ^ 2StatC^ for identical parti¬ 

cles (jOrmel fc Cuzzil[2007L note Birnstiel et al. drop the 
numerical factor), the largest particle in fragmentation 
equilibrium has Stokes number 


St, 


^ Q, 

2atc^ atc^ ’ 


(56) 


and c is nearly constant across the outer disk at least. 
Even in the inner disk where c varies somewhat more 
significantly, St is nearly constant, but at a smaller value 
commensurate with the lower Q, for silicates. 

In the Epstein flow regime, which applies to most of 
the particles sh own in Fig. El the Sto kes number can be 
written as fe.g.. IBirnstiel et d1l2011h 


St = 


S 


(57) 


For a roughly constant Stokes number as a function of 
semimajor axis, r oc E implying that the decrease in r, 
and tm in the outer regions simply mirrors the decrease 
in gas surface density. However, the behavior in the in¬ 
ner nebula is not reliably predicted from these formu¬ 
lae because the particle sizes are not strictly in the Ep¬ 
stein regime, but often in the transitional or full Stokes 
flow regime (which is generally ignored in other models). 
Though we cannot easily find simple expressions for ex¬ 
pressing St in terms of the transitional formula we use, 
we can still give estimates for St when r > Amfp (see 
Appendix A.4): 


16 Pp rc 

CdAEpg “ 

Rep < 1; (58) 

\ 1/2 

where in the case Rep > 1 the particle-to-gas relative 
velocities AFpg are assumed to be due to turbulence. In 
this work, the largest particles we grow have near unity 
particle Reynolds numbers. In general, for relatively 
smaller particles, the assumption of turbulent driven rel¬ 
ative velocities is probably adequate unless at is very 
small. However, as particles grow larger or in the outer 
nebula, their relative velocities may be driven by head¬ 
wind effects and closer to ~ rjVK- Later in Sec. Sjwe will 
compare these estimates to our simulations directly. 

Also in the middle panel of Fig. ID] the particle inter¬ 
nal densities are plotted for both the migrator and dust 
particles as a function of semimajor axis. The “cosmic 


2 Ppr^VK 
9 

1 

4 (ay^E)i/3 


abundance” particle density of ^ 1.5 g cm“^ is main¬ 
tained in the outermost nebula, but significant changes 
in mean particle density have occurred in regions where 
EFs are or were found. The smoothness of the particle 
density profiles reflects the process of sublimation from 
drifting particles at EFs, followed by outward diffusion, 
recondensation, and advection. The most striking ex¬ 
ample of this process is evident in the region between 
5 — 9 AU where the density of local solid material has 
essentially dropped to the density of water ice. The ra¬ 
dial extent of this water-ice enhanced region is partly a 
result of outward diffusion of vapor and advection of icy 
grains, but (in this case) also of the inward evolution over 
time of the water ice EF to its current location at ~ 4 
AU. The magnitude of this enhancement may be slightly 
overestimated if some amount of the drifting water ice 
were not evaporated inside the EF, but somehow buried 
and trapped in migrators to be carried further inwards. 
Finally, we note that slight differences in the density of 
migrators versus the “dust” population outside of EFs 
are due to the fact that we assume that condensation of 
vapor follows the surface area of particles (dominated by 
the smaller particles). 

For the same case “F”, Figure j^shows the evolution of 
gas (black curves) and solids (red curves) surface density 
at three different times (top panel), and (at 2 x 10® years) 
the radial velocities of solids and gas (middle panel) and 
the gas and solids accretion rates (lower panel). The 
gas radial velocities (middle panel, solid curve) are rel¬ 
atively low due to the choice of at, and change from 
inward {Vg < 0) to outward flow {Vg > 0) at around 7 
AU (see below). Also in the middle panel we plot the ra¬ 
dial velocities of the largest particle Vd,L (dotted curve, 
corresponding to the dotted curve in the middle panel of 
Fig. |3|) which shows how quickly in most places these 
particles are drifting inwards. Variations in Vd,L are due 
to differences in the size (and mass) of particles as well 
as transitions through EFs (Sec. I2.4.5|) . We also plot the 
radial velocities of the inward (V)j“, long dash) and out¬ 
ward (V^j*", short dash) drifting components of the dust 
population (Sec. 12 .2.21) . The (small) outward flow of the 
smallest grains (due to advection and diffusion) explains 
why there are solids further out than the initial 100 AU 
solids cutoff. The inward drift is strong in the outer disk 
and explains how the outer disk has been cleared of most 
of its material. We note that decreases sharply in the 
inner disk and eventually becomes similar to Vg. This de¬ 
crease can be associated with the drop in St inside the 
water ice EF, and leads to the subsequent pileup of solid 
material in the inner nebula. 

Reversal of the radial mass flow driven by viscosity oc¬ 
curs around 7 AU, with outward flows beyond this value, 
but in spite of this, by 2 x 10® years the bulk of the 
outer nebula solids have collapsed to within this radius. 
This can be seen in the surface density plot (top panel) 
where although the gas surface density has not evolved 
considerably due to the low at, the changes in the solids 
over time are pronounced (it is even more apparent in 
the lower panel of Fig. Ej). At 10^ years there is still 
a considerable amount of solid material out to 100 AU, 
but by 2 X 10® years, peaks in Ed outside the organics 
and water ice EFs are well defined, supplied by inwardly 
drifting material. Inside 7 AU, the advection of material 
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is inward, and although some vapor may diffuse outward 
to recondense, this appears to remain well enough inside 
the radius of mass flow reversal so that there is no longer 
an effective outward transport mechanism and material 
will remain confined to the inner nebula. 

Other effects are more subtle. In the lower panel of 
Fig. m sharp decreases in radial drift and an ensuing 
pileup of solids at the water ice EF are associated with a 
spike in the total accretion rate of solids Md (red curve). 
Just inside the water ice peak one finds that there is 
a bump in the gas accretion rate that corresponds to a 
local minimum in the pressure gradient (dashed curve). 
The radial drift of the lucky particles drops at 1.5 AU; 
this is associated with a local decrease in rn, related to a 
peak in the total migrator mass and the locally increased 
erosion that results. 

3.3. Fragmentation and bouncing (“BF”) 

When the so-called “bouncing barrier” is included in 
our fiducial model, the situation changes considerably. 
Because the particles are growing more slowly and re¬ 
main smaller (see below) the opacity is larger and the 
temperature in the inner nebula remains high even af¬ 
ter 2 X 10® years (see Figure [5l top panel), remaining 
more or less unchanged from its value at 10® years. Also, 
the water ice EF is found further out, at around 7 AU. 
The opacities are still decreasing with time in the outer 
nebula and increasing in the inner nebula, but remain 
generally larger than in the pure fragmentation case F. 
The fractional masses (lower panel. Fig. [5]) again show 
enhancement peaks in the solids near the organic and 
water ice EFs corresponding to local peaks in the opac¬ 
ity, but not so much the silicate EF. Between 10® and 
2 X 10® years, the EFs evolve only slightly in this case, so 
unlike the fragmentation-only case F, the enhancements 
in the solids are higher outside EFs - in fact, exceeding 
the enhancements of the vapor for the organics EF. The 
total enhancement of condensibles, in both solid and va¬ 
por form - in the inner regions is as high as ^ 5 — 6 
(black curve) by 2 x 10® years. The overall fractional 
mass of migrators is similar (~ 10“®) to the fragmenta¬ 
tion case, but they are now present everywhere in the 
innermost regions. A sharp cutoff in the migrator popu¬ 
lation around 5 AU is because particle sizes have not yet 
reached the fragmentation barrier outside this location. 
Interestingly, though no peak in the dust outside the sil¬ 
icates EF is quite discernable, the migrators do show a 
noticeable enhancement. 

Particle sizes plotted in the middle panel of Fig. [S] 
demonstrate that in most of the disk, the largest particle 
has not reached the fragmentation size (rn < r*) - so the 
heavy dashed curve for r* does not appear. The bouncing 
barrier has slowed the growth rate considerably because 
sticking becomes zero for equal mass particles at much 
smaller sizes than in the fragmentation-only case. In gen¬ 
eral, the bouncing barrier (BB) restricts eligible collision 
partners for growth to a narrower range of smaller pro¬ 
jectile sizes - causing all particles to grow more slowly. 
Still, as long as a population of small particles remains, 
the BB is not impermeable and we expect that the frag¬ 
mentation barrier will be reached at later times. 

After 2 x 10® years, only particles inside the water ice 
EF have reached r* (heavy dashed curve appears) which 


is understandable since the larger Q* outside this EF 
means it will take much longer to reach the larger frag¬ 
mentation limit r* (c/. Fig. [11] Sec. 13.51) . It is apparent 
that r*, tm , and rn in the inner nebula are smaller in this 
case than in case F. This seems counterintuitive since the 
strength remains the same; however, regarding tm and 
tl, the reason is that the number density of dust parti¬ 
cles is much higher and the probability of destruction for 
particles of size tm and rn is proportional to the number 
densities of the smaller particles they are colliding with 
(see Eq. [53]). StM is lower than St, which we argue 
is due to a slow growth or drift related effect which we 
discuss in Sec. KT[ 

The smaller and inwardly decreasing value of r* is ex¬ 
plained in a different way. The Stokes numbers (blue 
curves) are smaller by an order of magnitude than in 
case F. A smaller St* is explained by the higher opacity, 
temperature, and stronger radial temperature gradient 
than the lower and constant temperature in case F at 
2 X 10® years (see Eq. [56], also see discussion. Sec. II, 
and this leads to a smaller r,. The bulk of particles in 
the inner nebula have StM ^ 10“® causing them to drift 
at a much slower rate, and to be retained for longer pe¬ 
riods before being lost with an ensuing increase in mass 
enhancement. Outside the water ice EF, we see a similar 
trend in StM as in case F, but at smaller values because 
the fragmentation size has not been reached. 

Finally, the internal density of dust and migrator par¬ 
ticles (red and magenta curves, middle panel) retains in¬ 
teresting variability. In general, the internal density pro¬ 
files are characterized by decreases just outside the var¬ 
ious EFs; this is due to enhancement in the evaporating 
species as it diffuses back across the EF and condenses 
on to dust and migrator particles. The radial range over 
which the density decreases is due not only to vapor con¬ 
densation but to the diffusion and advection of smaller 
dust grains. The interplay is quite complex though, as 
we will see below from gas and dust accretion rates. In 
particular, the density variation outside the water ice EF 
for case BF does not span the same radial extent as in 
case F, because the EF has been relatively stationary 
over the simulation time in case BF. However, its edges 
are less sharp because of the more significant outward 
advection and diffusion of smaller icy grains. 

The radial velocities of the various particle sizes plot¬ 
ted in Figure [6] (middle panel) are lower than for case 
F, reflecting slower radial migration of smaller particles. 
The location where V)j“ becomes similar to the gas radial 
velocity is slightly further outward due to the smaller 
particle size, and in fact occurs just outside the organics 
EF. The inward drift of lucky particles reaches a min¬ 
imum speed at this point, and increases inwards up to 
the silicate EF which is at the inner edge of the compu¬ 
tational grid. The more complex behavior seen in the 
dust and gas accretion rates compared to case F again 
emphasizes that there is inward advection of material in 
the inner regions, although at speeds less than the gas 
velocity, but also outward diffusion at EFs further af¬ 
fecting the accretion rate. The gas accretion rate shows 
variation through the organics and water ice EFs, be¬ 
cause particles containing most of the mass are small and 
traveling with the gas. The mass accretion rates of both 
the gas and dust become more similar at these points 
implying that the dust is apparently influencing the gas 
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BF: vor Q„ = 4x10"^ 



Figure 5. Fiducial model with bouncing and fragmentation (BF). 
Panel descriptions same as Fig. 0 When the bouncing barrier is 
included in our simulations, particles grow more slowly and remain 
relatively small over the period of the simulation (r > r* inside 5 
AU only). Lucky particles are also smaller compared to the F case. 
As a result the temperature remains hot in the inner nebula, and 
opacities are generally larger there as well. Particle growth in the 
outer disk is still sufficient to lead to significant inward transport, 
though much less than the F case. Smaller particles also lead to 
relatively large enhancements in dust and vapor in the inner nebula. 
See text for details. 

evolution even though the solids-to-gas-ratio is not near 
unity. The reason is probably because opacity peaks lead 
to temperature gradient fluctuations, that lead to radial 
viscosity variations through the sound speed c, and thus 
corresponding variations in the gas radial flow. 


3.4. Fraqmentation, bouncinq and mass transfer 
(“MTBF”) 

The process of mass transfer (in high velocity col¬ 
lisions) encompasses an abbreviated range of condi¬ 
tions (Sec. I2.4.4L but can have a significant effect 
on certain aspects of particle growth over long periods 
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Figure 6. Fiducial model with bouncing and fragmentation (BF). 
Panel descriptions same as Fig. [4| Because growth times are 
longer, and particles generally smaller, the dust surface density 
has evolved much less than the F case after 2 x 10^ years. The ra¬ 
dial velocities of the particles sizes of interest all reflect this trend. 
The dust mass accretion rate is larger because more mass is in the 
dust, and shows substantial variation due to interaction with the 
various EFs. See text for details. 

(|Windmark et al.ll2012allbl . l2013t[Garaud et 3111201311 . In 
Figures 0 and [8] we show the results of our fiducial 
model that also includes mass transfer. The evolution 
of temperature (top panel, Fig. somewhat differs be¬ 
tween this case and the previous case BF, in that T in 
the innermost nebula is slightly higher by 2 x 10® years, 
even though the water ice EF has moved inwards. The 
opacities also are similar in magnitude, but their radial 
distribution is different. The water ice EF has a larger 
radial extent here. Compared to the BF case, there is a 
much larger drop in the opacity inside the ice EF and a 
much larger opacity enhancement outside it, suggesting 
(like in case F) that there is a larger difference between 
particle sizes inside vs. outside the water ice EF. The 
large contrast across the water ice EF is also visible in 
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MTBF: vor Q„ = 4x10"* 



Figure 7. Fiducial model with mass transfer, bouncing and frag¬ 
mentation (MTBF). Panel descriptions same as Fig. [S] In general, 
the MTBF case is quite similar to BF with the exception that 
growth in the migrators is slightly faster, and particle sizes are a 
bit larger. Overall, the enhancements one finds in both the dust 
and vapor in the inner nebula are similar to the BF case, except 
that more mass appears in the dust and less in the vapor. The 
temperatures in the silicate rich region are slightly hotter as well. 
See text for details. 

the fractional mass plot (lower panel). The outer disk 
remains similar to case BF, but differences in the inner 
regions can be associated with the evolution of the par¬ 
ticle size distributions. Overall, the total enhancement 
in the inner nebula is similar to the BF case, a factor of 
^5 — 6, though with slightly more in the solids compo¬ 
nent and less in the vapor. 

The contrast between particle sizes inside and outside 
the water ice EF is shown in the middle panel, where 
compared to case BF, growth at the larger sizes is faster 
and as a result the fragmentation size r, has been reached 
further out in semimajor axis by 2 x 10^ years. The ef¬ 
fect of mass transfer appears to lead to two things. A 
slightly larger value of tm relative to r* most easily seen 
inside i? < 1.2 AU, and consistently larger lucky parti- 
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Figure 8. Fiducial model with mass transfer, bouncing and frag¬ 
mentation (MTBF). Panel descriptions same as Fig. [4] The overall 
trends are similar to the BF case, especially in the outer nebula. 
Differences in the inner nebula can be attributed to the evolution 
of the particle size distribution when mass transfer is included. See 
text for details. 

cles. These effects are much more evident when compar¬ 
ing the middle and lower panels of Fig. [5] which shows 
particle size distributions over this radial range. The 
lower panel, which refers to the MTBF case, is character¬ 
ized by secondary peaks which are indicative of particles 
that are benefitting from particle collisions in which the 
growing particle accretes more mass than is eroded (see 
Sec. 12.4.41) . A slight dip in rp appears outside the the 
water ice EF which is associated with the large spike in 
fractional mass there. The particle mean density profiles 
are not much different from the BF case, except in the 
water ice enhanced region, with the differences tied to 
the different radial extents of the organic and water EFs 
in the two cases. In the ice-enhanced region just outside 
the ice EF, both tm and rp are larger when mass transfer 
is playing a role. 

The velocity trends in Figure [8] are qualitatively sim- 
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Q. = 10^ Oi = 4x10-'* 



Figure 9. MTBF model assuming a constant fragmentation en¬ 
ergy Q* = 10® (Q103). Panel descriptions same as Fig. |3] 

ilar to the BF case. The most noticeable differences ap¬ 
pear in the radial velocities (middle panel). The pres¬ 
ence of much larger lucky particles inside the water ice 
EF that are migrating inwards with relatively larger ra¬ 
dial velocities (Eq. [12]), accounts for the much steeper 
inwards increase in Fd,L compared to the BF case. Also, 
whereas drops off gradually to its lowest value in the 
BF case between ~ 3 — 4 AU, a relatively sharp bump 
appears in the MTBF case because the fragmentation 
barrier has been reached and most of the mass here is 
in larger particles which influence the gas. This coin¬ 
cides with a stronger peak in M than what was seen at 
the same location in the BF case. The corresponding 
pressure gradient also a takes a much sharper drop here 
(lower panel). 

3.5. Models with constant Q* 

We have employed a variable Q* (and bouncing bar¬ 
rier) in all of the models presented thus far, motivated by 
experimental results that suggest that icy aggregates are 
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Figure 10. MTBF model assuming a constant fragmentation en¬ 
ergy Q* = 10^ (Q103). Panel descriptions same as Fig. [4| 

stickier and indeed stronger, at the velocities in qnestion, 
than silicate aggregates (see Sec. 12.4.21) . In this section, 
we present MTBF simulations that assume constant val¬ 
ues for material strength, choosing two extreme values 
for Q* = 10^ and Q* = 10^ which would correspond to 
fragmentation velocities of ^ 30 cm s ~^ and ^ 30 m s~^ , 
respectively. Of these two, results of iBeitz et al.l (|201lL 
for silicates) suggest that the larger value is especially 
unrealistic. _ 

In Figures [9] and 1101 we show our simulation re¬ 
sults for Q, = 10^. Not surprisingly, this small value of 
Q* means that particle sizes (middle panel. Fig. [9]) are 
much smaller throughout the entire disk. As an example, 
r* ~ 0.18 cm at 1 AU in the MTBF case (where Q, = 10^ 
for silicate grains) whereas r* ~ 0.02 cm for this case, 
nearly an order of magnitude smaller. The maximum 
fragmentation size occurs around 10 AU. Because parti¬ 
cles remain smaller, Stokes numbers are smaller which 
means there is less radial drift (and loss) of material 
over 2 x 10® years. The Stokes numbers for tm decrease 
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Q. = 10^ Oi = 4x10-'* 



Figure 11. MTBF model assuming a constant fragmentation en¬ 
ergy Q* = 10^ (Q107). Panel descriptions same as Fig. |3] 

monotonically from larger values in the outer disk to 
StM < 10“^ in the innermost regions, and we do not 
see the roughly constant values as we did in the previous 
cases. Interestingly, however, we do see the Stn roughly 
constant, or r oc S, relationship hold for the lucky parti¬ 
cles in the outer nebula. 

The temperature and opacity (upper panel) and the 
fractional masses (lower panel) further emphasize the 
sluggish disk evolution this case results in. The disk 
has only cooled significantly in the last 10® years, but 
still remains relatively hot partly because of a gradual 
increase in inner nebula opacity due to a small enhance¬ 
ment of solids. Unlike other cases we have seen previ¬ 
ously, the outer edge of the water ice EF has moved very 
little, with a slight increase in radial extent. Because 
the fragmentation barrier is at such small sizes, migra¬ 
tors appear throughout the disk at a fraction of ^ 10% 
of the dust component. Moreover, the fractional mass 
in the outer disk has remained similar to its initial value 
because the gas and dust are still well coupled. Despite 
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Figure 12. MTBF model assuming a constant fragmentation en¬ 
ergy Q* = 10^ (Q107). Panel descriptions same as Fig. 0 

the sharp pressure gradient at the outer edge (dotted 
curve, Fig. 1101 lower panel) a significant amount of solid 
material has advected outward, beyond 100 AU. The ra¬ 
dial velocities in Fig. [10] (middle panel) are very close to 
the gas velocity everywhere except approaching the disk 
edges. In fact the inward drifting dust velocity compo¬ 
nent Vy diverges from nearly zero only at the outer edge 
of the nebula solids disk. The surface densities (upper 
panel) also show generally little evolution, and the dust 
accretion rate (lower panel) is muted. A low Q, model 
thus appears to be a way to minimize the amount of 
sol ids loss over long tim escales, as previously suggested 
bv IBirnstiel et al.l (120091) . 

Figures 1111 and 1121 for Q* = 10^ show an entirely 
different behavior. In this case, the fragmentation size 
is so large that even after 2 x 10® years, r < r, ev¬ 
erywhere (heavy dashed curve does not appear), even 
though tm is as large as ~ 10 cm in the inner neb¬ 
ula (Fig. [TTl middle panel, dotted curve). The Stokes 
number neatly follows the nearly constant theoretical 
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LI: vor Q., = 4x10-'* L1: var Q., = 4x10-'* 



Figure 13. MTBF model assuming a star with 1 solar luminosity 
(LI). Panel descriptions same as Fig. [3l 

“equilibrium” value in the outer disk, giving r oc E (see 
Sec. 13.2|) but diverges from this relationship in the in¬ 
ner regions possibly because these relatively large parti¬ 
cles may be in the transitional Stokes regime. The disk 
temperature (top panel) decreases quickly because the 
particles grow relatively quickly, and much of the ma¬ 
terial in the inner disk is beginning to be lost, leading 
to the substantial decrease in opacity. Despite the sys¬ 
tematic loss of material, we still see two strong peaks 
in opacity due to the enhancements outside the organ¬ 
ics and water ice EFs. These are also mirrored in the 
fractional masses (lower panel). We note an apparent 
peculiar difference between this case and other cases we 
have shown so far: the outer nebula solids disk has not 
“collapsed” as far as other models, despite the fact that 
the higher Q* here allows growth to much larger parti¬ 
cles, which should drift more quickly. The explanation is 
that the bouncing barrier we assume here (as in the case 
of Q* = 10^) is also a constant set at the silicate value, 
but is not scaled with Q* (see Sec. I2.4.2|) . so that the dis- 


Figure 14. MTBF model assuming a star with 1 solar luminosity 
(LI). Panel descriptions same as Fig. [S] 

parity between the bouncing and fragmentation barriers 
is especially considerable in this case. Not only does this 
mean that growth is slower since larger particles can only 
grow from particles smaller than the bouncing size, but 
the number density of those bouncing barrier-or-smaller- 
sized particles is always decreasing here because they are 
not being replenished by fragmentation. 

Finally, looking at the radial velocities in Fig. [12] (mid¬ 
dle panel), we see that the radial drift rates are larger in 
the inner nebula for both the largest particles (Vd^n) and 
the dust (V)j“), explaining the significant loss of material 
from the inner nebula, the opposite of what was found 
in other cases. In this case, the gas surface density (top 
panel) has evolved the least of all the cases we have shown 
thus far, with values very close to the initial profile even 
after 2 x 10^ years. In fact, in the innermost region, 
the gas density profile has become even steeper. This is 
because the rapid growth and loss of inner solar system 
solids has led to the lowest opacity, producing the fastest 
cooling and lowest viscosity. Remember that results like 
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this depend on the Ot-model prescription for viscosity 
increasing with temperature. The steady inward migra¬ 
tion of material is also evident from the dust accretion 
rate that remains positive throughout the disk. The gas 
accretion rate is much less strong in this case in the inner 
regions because of the rapid cooling and steeper surface 
density. 

3.6. A 1 Lq model 

We have explored an MTBF model in which we fix the 
luminosity at 1 Lq , as presented in Figures 1131 and 1141 
Previous models have all used the time dependent lumi- 
nosity model of iD’Antona fc Mazzitellil (Table 3, 11994) 
where the initial luminosity is of order ^12 Lq and de¬ 
creases to ~ 3 Lq after 1 Myr. Thus here we can compare 
the effects of luminosity in the MTBF model. The pri¬ 
mary effect one would expect is on the disk temperature 
(top panel of Fig. [T31). We find only subtle differences 
in the overall temperature profiles, which are primarily 
associated with the radial extent of the water ice and or¬ 
ganics EFs, and the overall steepness of the temperature 
profile inside this location. In particular, the gradient 
is so steep that the organics EF is traversed within one 
radial bin. The inner nebula shows only small changes 
because its energetics is dominated through these stages 
by local viscous dissipation rather than stellar irradia¬ 
tion. Differences in the Rosseland mean opacity profiles 
suggest that growth is proceeding a bit more rapidly in 
the outermost portions of the disk, at least at the earliest 
stages (10^ years). This is also evident in the fractional 
mass plots (lower panel). The 10® year fractional mass 
curves look quite similar for this case and the fiducial 
MTBF case, but have slightly larger radial extent be¬ 
cause of the location of the EFs. However, by 2 x 10® 
years the vapor and dust content in the inner regions 
are considerably higher than in the MTBF case, with an 
overall enhancement of an order of magnitude over the 
initial profile. The Stokes numbers corresponding to tm 
and tl (middle panel) are even flatter across the outer 
disk (and higher than in the MTBF case, which may 
explain why the outer disk evolved faster in the early 
stages here). In the inner nebula, there are also differ¬ 
ences in the particle sizes relative to the MTBF case, 
which may derive from the increased number density of 
particles making mass transfer only slightly more effec¬ 
tive for particles growing from r* to near the size of cm. 

In Figure 1141 the gas accretion rate profile is flat¬ 
ter, and its magnitude lower than the MTBF case. This 
is once again due to the steeper temperature gradient. 
In particular, one finds that the pressure gradient plum¬ 
mets much more sharply through the organics EF with 
this location corresponding to the sharp decrease in the 
solids accretion rate at ^ 2 AU just inside the strong 
organics peak. The peak in Md near the water ice EF 
is more intense, meaning that more solids are being de¬ 
livered to the inner disk at this time compared to the 
MTBF model. Although a higher StM in the outermost 
part of the nebula implies faster radial drift rates, these 
are hardly noticeable in the radial velocity plot (middle 
panel). In general, the radial velocities are similar to the 
MTBF case. Overall, the 1 Lq model makes little differ¬ 
ence to disk temperature evolution at these early stages, 
though the outer disk is colder, and this may have led to 


the subtle changes that we see in this model compared 
with the variable (higher) luminosity MTBF case. At 
later times, the higher (more realistic) luminosities may 
play a more important role in a relative sense. 

3.7. MTBF model with high at (“HA-MTBF”) 

As a final case, we compare a more viscous disk to 
our fiducial MTBF model; here the turbulence parame¬ 
ter at = 4 X 10“® iFigures 1151 and I16p . We expect that 
regions which are dominated by viscous heating should 
be hotter in this model, and indeed Fig. [15] (top panel) 
shows high temperatures extending to larger distances 
from the star at the nominal time of 2 x 10® years. It 
is interesting that the higher at disk starts hotter than 
the nominal MTBF disk and cools to the current values, 
whereas in the fiducial MTBF model (and LI as well) 
the temperature cools and then actually increases to the 
value at 2 X 10® years. The greatest differences we see 
between the two values of at, apart from the full silicates 
EF being present, have to do with the radial extent of 
both the organics and water ice EFs. Both fronts are fur¬ 
ther out, but are also smaller in radial extent than the 
MTBF model. The systematically decreasing T is natu¬ 
rally due to the decrease in Rosseland mean opacity as 
particles grow. Yet, the MTBF case shows an increasing 
T despite having larger particles because the enhance¬ 
ment leads to a higher opacity. In the HA-MTBF model, 
there is very little enhancement as is seen in the fractional 
mass (lower panel). Another clear effect of the more vis¬ 
cous disk is that the nebula gas evolves more strongly 
outwards, carrying a significant amount of solids with it 
(Fig. m upper panel). The strong outflow of solid ma¬ 
terial is also evident in Fig. [T5| (lower panel) where after 
2 X 10® years, the edge of the solids disk extends to about 
^ 300 AU. Although the fraction of migrators is compa¬ 
rable between this and the MTBF case, the dust fraction 
in the innermost disk is less by an order of magnitude. 
Most of the solids are in the vapor phase. 

The particle sizes (middle panel) also differ in an un¬ 
derstandable way from the fiducial MTBF model. The 
relative velocities between particles are larger for the 
larger at, which means that the bouncing and fragmenta¬ 
tion barriers are both reached for smaller particles. Com¬ 
paring the middle panel of Fig. [13] with that of Fig. [T] 
we see that the fragmentation size r, is now over an or¬ 
der of magnitude smaller. Interestingly, even while the 
mass-containing particles are smaller here, the lueky par¬ 
ticles are not too much different in size from those in the 
hducial MTBF case - so particles here can get even luck¬ 
ier relative to the nominal fragmentation barrier! This 
again may be due to the radial drift barrier being the 
limiting factor in how large these particles can grow (see 
discussion on Stokes numbers. Sec. Ej). Increasing at by 
an order of magnitude has led to a smaller StM in the 
inner regions - also by an order of magnitude, consistent 
with Eq. (EH). However, StM values are hardly changed 
outside the water ice EF. This suggests that the inner 
nebula is obeying fragmentation-barrier limitations, but 
something else - probably radial drift - is playing a role 
in the outer nebula (see Sec. 14.11) . 

In Figure [m top panel, we can see that the larger at 
for this case has evolved the gas density much more, with 
E already at values typical of or sm aller than a minimum 
mass nebula (e.g., iHavashil 119811) outside ^4 — 5 AU, 
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HA-MTBF: vor Q,, a^ = 4x10"^ 



Figure 15. MTBF model assuming at = 4 X 10 ® (HA-MTBF). 
Panel descriptions same as Fig. |4] 

while inside the evolution has not proceeded as quickly 
and is similar to the MTBF value at 2 x 10^ years be¬ 
cause of the hotter, steep temperature gradient disk. In 
fact, the accretion rate of the gas (lower panel) shows 
that M lingers close to zero between ~ 1 — 4 AU indi¬ 
cating very small gas advection velocitie^. The mass 
flux reversal occurs much further out in the disk than 
the fiducial MTBF model, as one might expect with a 
much more viscous disk. The variations in are due 
to the various EF transitions but are more muted than 
the fiducial MTBF case again because of a low level of 
mass and opacity enhancement. The ‘bumps’ in solids as 
seen in the fractional mass profile may still be building 
up, however. Finally, the radial velocities (middle panel) 

This is mainly due to the steepness of the temperature gradient. 
From inspection of Eq. we can see how can become positive. 
For example, if the temperature T oc R~^, it is easy to show that a 
reversal in sign will occur if, ignoring a term in the surface density 
oc dlnU/dlnR < 0 for simplicity, 2 — s < 0, where s is the power 
law exponent. 
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Figure 16. MTBF model assuming at = 4 X 10 ® (HA-MTBF). 
Panel descriptions same as Fig. |4] 

demonstrate that outward transport is more significant 
than the lower viscosity models. The gas radial velocity 
Vg increases more steeply in the outer nebula, and this 
affects the sizes (and amount) of particles that have pos¬ 
itive (outward) velocities I/j"*". The inward drifting solids 
component Vy is stronger than in the fiducial MTBF 
case, suggesting that at later times, vast regions of the 
outer nebula (between the water ice EF and ^ 100 AU) 
the disk may become depleted in solids. 

4. DISCUSSION 

In this section we review the effects of different stick¬ 
ing/bouncing/fragmentation schemes on particle size dis¬ 
tributions, and note some implications regarding mm-cm 
wavelength o bser vations of extrasolar disks (Sec. 14.11) . 
Also in Sec. 14.11 we digress into how the tradeoff be¬ 
tween growth and drift in evolving nebulae may explain 
a small but interesting difference between our results and 
previous work. Next, we look more closely at general im¬ 
plications of growth and drift for radial redistribution of 
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Figure 17. Summary of a the particle size distributions after 2 X 10® years for a suite of models we have presented in this paper. 


mass (Sec. 14.21) . and note implications for meteorites, 
planetesimal formation, and perhaps even the radial dis¬ 
tribution of planets. In Sec. 14.31 we look more closely at 
the effects of Evaporation Fronts, specifically on produc¬ 
ing large objects enriched in, e.g., water ice. We next 
discuss implications of our self-consistent opacities, in 
particular for mechanisms of turbulence generation and 
nebula gas evolution (Sec. 14.411 . We conclude with some 
caveats and a short discussion of future work. 


4.1. Particle size distributions; the drift and bouncing 

barriers 

In Figure 1171 we collect the particle size distribution 
results after 2 x 10® years for the suite of models described 
in Sec. |3l Though there is considerable variation in 
particle sizes across the models, there are specific trends 
that are characteristic of all of these results. 

In all figures, particle radii are shown in black and the 
corresponding Stokes numbers in blue. Particle internal 
densities are shown by the red/magenta curves. Recall 
that r* (dashed lines) is the fragmentation radius, and 
rM (solid lines) is the radius dominating the total particle 
mass. When r* and tm are plotted, the dotted line (rp) 
refers to the largest “lucky” particle in a radial bin. Since 
particles of radius rp contain a negligible amount of mass 
fsection |3.I.lL we emphasize the mass-dominant particle 


rM- Recall that, if the heavy black lines for r* and ryi 
do not appear, it indicates that the largest particle tl 
( which in this case does contain most of the mass) has 
not yet reached the fragmentation size r,. 

The critical role of composition-dependent strength Q, 
for limiting particle size is clear. In all the nominal cases 
where we allow icy particles to be stronger than ice-free 
particles, a matching step-function is seen in tm and StM- 
The constant Q, cases, naturally, do not show this step. 
This simple result is explained, to first order, by the sig¬ 
nificant change in the fragmentation barrier (equation 
[56l) when Q* changes so significantly (see below for more 
discussion). In section 02] we discuss how this size dif¬ 
ference leads to strong mass enhancement /depletion dif¬ 
ferences between the inner and outer nebula. Another 
general result is that while “lucky” particles tl ^ tm 
can be found throughout the inner nebula (although car¬ 
rying negligible mass), the outer nebula is devoid of such 
particles. That is, in the outer nebula the size distribu¬ 
tion cuts off more sharply at its largest size. The same 
result is manifested in the Stokes number plots; this is 
discussed below. The most realistic cases (MTBF and 
HA-MTBF) show a peak in StM and tm just outside the 
ice EF, probably because of the locally enhanced solids 
mass. However, this does not hold true in the inner neb¬ 
ula where there is little if any correlation with the inner 
EFs and the distributions in tm and StM, with perhaps 
the exception of case F. The near-constant values of StM 
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in most of the plots in Fig. [T7] suggest fragmentation- 
dominated regimes, as can be understood from equation 
[56] in sectioning (recalling that T is only weakly varying 
in, at least, the outer nebula). 

For our most realistic models (BF and MTBF), the 
inner nebula where the asteroids and terrestrial plan¬ 
ets form (where the particle strength is lower because 
ice is absent) is capable of supporting growth only to 
0.1—several mm r adius, as has been shown previously 
(jZsom et al.ll20inf l. This size is typical of “chondrules” 
which dominate unmelted meteorites. Our particles are 
not themselves “chondrules”, which most workers believe 
have been melted as freely-floating nebula objects, but 
represent aggregates of the right size to serve as their 
precursors. Just outside the ice EF, radii range from 
several dm to almost a meter. This might seem sur¬ 
prising in view of the iconic “meter-size barrier”, which 
is generally placed at St ^ 1 in a Minimu m Mass So¬ 
lar Nebula (iCuzzi fc Weidenschillind [200611 . but notice 
from the Stokes number curves that StM remains in the 
range 0.01 — 0.1 because at this early stage the nebula 
gas density remains fairly high. This means that radial 
drift velocities and collisional velocities are still increas¬ 
ing as these particles grow. The sizes of particles in the 
outermost nebula depend most strongly on the strength 
of turbulence and the assumptions of the sticking model. 
Around 30 AU for case F, particles can reach 1 cm radius, 
but for other, more plausible cases, sizes at 30 AU are 
closer to 1 mm because of the inhibiting role of bouncing. 
Thus, even at this early stage, “bouncing” has played an 
important role in limiting the sizes of particles detectable 
by mm-submm observations (for more discussion see “Di¬ 
gression” below). 

Our results thus ge nerally agree with the results of 
previous workers (e.g., iBrauer et al.ll2005lBirnstiel et al.l 

120091 . [2?)Tol : lZsom et al.ll20inl) that the fragmentation and 

bouncing barriers (along with radial drift; see below) 
limit particle growth. When only fragmentation is in¬ 
cluded (F), growth is quite rapid and stalls at particle 
sizes that are the largest amongst all of these cases. Yet, 
because particle growth reaches larger sizes, the disk 
mass drains most rapidly into the inner nebula due to 
radial drift (see Fig. [31 lower panel), with important 
consequences, as discussed below. The bouncing barrier 
(BF case), even if mass transfer is added (MTBF case), 
effectively reduces the rate of growth because a target 
particle is limited in the size of particles it can grow from. 
Thus in Fig. [T3 we see that the fragmentation barrier 
has not been reached in the outer nebula by 2 x 10® 
years for the BF and MTBF case s. Like previous work¬ 
ers (e.g., iWindmark et al.ll2012alibl : iGaraud et al.l 1201^ 
iDrazkowska et al.l 120131 12014( 1 we find that incorporat¬ 
ing high-speed mass transfer (see Fig. [2]) can allow the 
growth of a small number of “lucky particles” to sizes 
tl > r,; however, the mass contained in this population 
is extremely small at best, and they are absent in regions 
dominated by drift (see below). 

Figure 1181 compares several important metrics for a 
subset of these models to the observed, mass-dominant 
value StM (solid blue). As noted above, the near¬ 
constancy of the mass-dominant particle Stokes number 
StM through regions where Q* is constant is largely ex¬ 
plained by Eq. l|56l red solid curve for St*). However, 
while all the models follow the general trend of St*, and 
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Figure 18. Comparison of our estimates for the fragmentation, 
bouncing, and radial drift Stokes numbers to the computed Stokes 
numbers for models F, BF and Q103 after 2x 10^ years. Radial drift 
dominates the outer disks in models F and BF. St* > Stjvi because 
systematic velocities are comparable to the turbulent velocities for 
these sizes. Only in model Q103 is > tgr in the outer disk. See 
text for details. 

quantitative agreement is excellent for the moderately 
realistic F, BF, and MTBF cases in the inner nebula 
(and at all radii under conditions of weak particles for 
Q* = 10^)7 quantitative agreement is poor in the outer 
nebula once bouncing and Mass Transfer are introduced 
(in the sense that St^i < St*). This disagreement arises 
because Eq. (EH) assumes that the collision velocity is 
dominated by turbulent relative velocities, and it turns 
out that, at least for the models presented here, given 
the small values of StM, the large values of rj in the outer 
nebula drive systematic drift-and-headwind related ve¬ 
locities that exceed the turbulent velocities, so Eq. (1531) is 
an overestimate of the fragmentation barrier size. When 
bouncing is included, the discrepancy (vertical offset) in¬ 
creases since approach to the fragmentation barrier is 
now slower (particles larger than rb can now only grow 
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from the small fraction of particles r < rb), and StM 
is “quenched” at smaller, non-equilibrium values by the 
evolving radial drift barrier, as described below. In the 
low-strength Q* = 10^ case, the fragmentation barrier is 
reached quickly everywhere with StM St*. In the inner 
solar system, for cases F, BF, and MTBF, turbulent rel¬ 
ative velocities dominate so the nominal St, expression 
appears valid in all three cases. However, looking deeper, 
the situation actually presents a subtle puzzle, and to un¬ 
derstand this we now estimate “limiting” Stokes numbers 
for bouncing and drift. 

Digression: role of the “Drift Barrier” in setting tm 
and ri^: We can derive a limiting condition for bounc¬ 
ing by setting Shim = m' = ms) = 0 in equation (pOl) . 
leading to = Co/nib- As in deriving the fragmenta¬ 
tion limit, we set Vb equal to the typical collision veloc¬ 
ity in turbulence for some radius rb: AFpp = c\/Stbat, 
where the particle size can be related, approximately, to 
the corresponding S tb (using Eq. (l57ll : also see, e.g., 
iBirnstiel et al.l[20ilh . We thus get the bouncing barrier 
Stokes number Stb: 


Stb 


iratC^T,^ j ’ 


(59) 


for which, in typical conditions, the strongest dependence 
is on the gas surface density. 

While the importance of a “radial drift barrier” has 
been known since the early work of Weidenschilling (1988 
et seq ), it has only recently been quantifi ed and its in¬ 
fluenc e is still being assessed. For instance, iBrauer et al.l 
(|2008ll simply assumed that St = 1 represented the radial 
drift b arrier; this was echoed by IBirnstiel et al.l (|2009L 
I20I0H . IBirnstiel et al.l (|20I2all used a more realistic ap¬ 
proach based on finding the particle size for which growth 
an d drift timescal e s were equal. We follow the approach 
of IBirnstiel et al.l (|20I2ari and generalize it in several 
ways. 

We derive our expression for the radial drift limited 
Stokes number Std by comparing the radial drift time 
td = R/Ur ^ i?/??VKSt to the particle growth time tgr = 
m/m. As in Eq. (PI). m ^ TTr^pdAV^gS*, where S' < 1 
is some average sticking coefficient which we will take 
as < 0.5 for specifics. Approximating AV^g ^ c-\/cqSt 
as the turbulent relative veloci ty, whether for same- 
size or different-size particles (lOrmel fc Cuzzil 1200711 . 
and the near-midplane density to be pd = Ed/2/id = 
fT,y'St/at/2H if = fd here), we get t„ = 2/3fQ, es¬ 
sentially the same as IBirnstiel et al.l (|20I2al their Eq. 
13). Then, we simply define the particle radius rd as the 
“drift barrier”, such that td ^ tgr: 


, 3^2 /SE 

’’d > -■ 

8 PpV 


(60) 


Note that the “drift barrier” is not simply determined 
by a drift velocity, but also depends on the growth rate 
through the amount of material available locally and the 
typical sticking coefficient. We relate rd to a correspond¬ 
ing drift-limited Stokes number Std, in either the Epstein 


or Stokes regime, using Eqns. (EZD or (|5S|) : 

^ (9/4)Amfp; 

f-Y Rep<l; 

(f) ^ Rep < 800, 

(61) 

If particles have St > Std, they are drifting faster than 
they can grow. The first expression, valid for the Epstein 
regime , corresponds to equation (17) of IBirnstiel et al.l 
(|2C)12all ; however, our expression includes a factor of 
S' < 1 and retains the full value of ry, which is not al- 
ways well repr e sented by <? /V^ as in equation (17) of 
IBirnstiel et al.l (|2012ari . In both cases, a higher local 
mass density makes it easier for growth to dominate drift 
(increases Std), and a stronger radial pressure gradient 
makes drift more dominant (decreases Std). 

We can now compare the calculated Stokes numbers 
(in blue), for the mass-dominant (StM; solid) and largest 
(Stb; dotted) particles with our estimated limiting “bar¬ 
rier” values for fragmentation (St*; red solid), bouncing 
(Stb, dashed red) and radial drift (Std, solid and dotted 
magenta as related to tm and tl respectively). As an 
indication of where the Stokes-Epstein regime transition 
lies, we also plot our approxi mation Eo. (l58ll for the 
Stoke s regime (dot ted red, see ICnzzi fc Weidenschillind 
120061 and Sec. 13.21) . Our actual models (blue curves) use 
a smoother bridging expression, averaging over all parti¬ 
cle sizes; the discrepancy between the red and blue dotted 
curves indicates that the effects of Stokes drag actually 
extend over a broad range of nebula radii. The bouncing 
barrier is active almost everywhere (Stb < StM); as men¬ 
tioned earlier, Stb does not provide an impermeable bar¬ 
rier, but merely slows growth for larger particles, which 
in principle could still reach the fragmentation barrier r* 
given a long enough time (in the absence of drift). 

The subtle puzzle in Fig. [18] is that while 
IBirnstiel et al.l (I2012al lbll found, in general, that the ra¬ 
dial distribution of the size of the mass-dominant par¬ 
ticles in the outer nebula corresponded very well to the 
radial variation of Std when drift dominated growth (ac¬ 
tually they present the corresponding rd instead of Std), 
our results show StM to have the shape of St* rather than 
that of Std, even when drift dominates (St > Std). 

We believe the explanation lies in a non-equilibrium be¬ 
havior to which the nebula gas evolutionary models pre- 
se nted here lend t h emselv es, but perhaps those adopted 
by IBirnstiel et al.l (l2012al lbll do not. Figure [T9| shows 
the contents of Fig. [T8|for case F, at three different times. 
The “drift barrier” Std is larger than StM throughout 
the nebula at early times, so growth dominates every¬ 
where. The fragmentation barrier is reached at larger 
distances as time goes on because at early times / is 
nearly constant, and growth is more rapid at smaller dis¬ 
tances where H is smaller (see tgr above Eq. l60l) . By 
2 — 3 X 10'^ yrs, St* is reached nearly everywhere, while 
Std > StM- Meanwhile, growth to large values of StM 
leads to rapid inward mass loss by radial drift, which de¬ 
creases /, which in turn decreases Std, so by 5 x 10^ yrs, 
Std has fallen well below StM in the outer disk and by 
10® years even for i? > 5 — 7 AU, and growth has ceased. 
There is, however, nothing to decrease particle sizes be- 
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Figure 19. Comparison of the fragmentation, bouncing, and ra¬ 
dial drift Stokes numbers to the computed Stokes numbers as a 
function of time, in case F. Here the transition from the growth 
dominated to the radial drift dominated regime is clearly shown. In 
a relative short time, the outer nebula solids are effectively trans¬ 
ported to the inner nebula. See text for details. 

cause SIm < St* still. Essentially the fragmentation- 
limited size of earlier times becomes frozen in, a fossil of 
an earlier era, while mass outside the ice EF continues to 
be lost to the inner nebula. A similar crossover behavior 
is seen for Std and StM in the BF and MTBF cases, at 
different times. 

This behavior has interesting implications for mm- 
submm observations of disks. If it is true that the 
particle size distribution faithfully mirrors the instan¬ 
taneous, underlying sur f ace ma ss density, as in the 
models of IBirnstiel et al.l (I2012al lbll. then measurements 
of size through spectral index determinations should 
be consistent with in ferences about total mass. Yet, 
IBirnstiel et al.l (|2010bD actually come up with inconsis¬ 
tencies when comparing models with observations; they 
found that ad hoc reduction in the surface mass density is 
required relative to the model. While case F (Fig. [T9l) is 


an extreme and indeed implausible case (the dominance 
of drift means that almost no mass is left over to form 
ice giants), the general decoupling of particle size distri¬ 
bution from simultaneous surface mass density provides 
a cautionary note. Other models (such as the MTBF 
case) also show an apparently fragmentation-limited size 
distribution represented by St, nearly constant, but also 
show significant mass de pletion relative to cosm ic abun¬ 
dance. We believe that IBirnstiel et al.l (l2012allbl l found 
a different behavior than do we because their underlying 
gas nebula, while formally evolving, is extremely large, 
so evolves in a much more leisurely way than does ours. 
It suggests that choice of initial disk size is something 
worthy of careful thought. As noted in section [2Tl differ¬ 
ent studies have adopted very different characteristic ini¬ 
tial nebula lengthscales, and these differences may have 
unforeseen but interesting implications. For instance, it 
might be that if one does measure the solids surface mass 
density and particle size distribution of a disk simulta¬ 
neously, one may be able to infer something about the 
evolutionary history of the region, and this is an exciting 
possibility. 

Finally, the changing role of Std vis-a-vis StM helps us 
understand the variable abundance of “lucky” particles 
(the variable behavior of Stn vis-a-vis SIm) as seen across 
all the models of Figs. HH] and UHl in a simple way. In 
all cases, lucky particles vanish when Std < StM- This is 
because growing “lucky” particles relies on long periods 
of inefficiently sweeping up much smaller particles at a 
slow rate. Particles simply do not reside long enough in 
a region where drift dominates to exceed StM, so their 
growth must have occurred entirely in environments en¬ 
countered earlier, which are no longer present. As they 
drift inwards, they do grow, although too slowly to reach 
the fragmentation barrier. Only in a growth-dominated 
regime where drift is not a factor, can “lucky” particles 
Stn > StM appear. Inspection of Fig. [12] indicates that 
this rule holds across all three cases (as well as others 
not shown here). Lucky particles are common (although 
a negligible fraction of the total mass) in the inner nebula 
in all cases, because Std > StM (Fig. (121). We note, how¬ 
ever, that the frequency and abundance of lucky particles 
may be signi ficantly increased if w e implement porosity 
(Sec. I2.4.ip . lOkuziimi et al.l (I2012D find that if the par¬ 
ticle porosity can get extrememly high (volume density 
of ^ 10“®) then drift behavior can change qualitatively. 
We explore this in future work. 

4.2. Radial redistribution of mass, and some 
implications 

We hnd that enrichment of the inner nebula in solids 
and vapor is a general, but not universal, result which 
seems dependent on the different Q, associated with the 
inner outer nebula. As we have seen in our simulations, 
the higher fragmentation energy of icy particles allows 
them to grow large enough to drift radially in significant 
abundance over relatively short time scales, even with 
the bouncing barrier imposed. The lower fragmentation 
energy of silicate particles in the inner disk means that 
particles originating in, or drifting into, that region, af¬ 
ter losing their ice, can only reach or maintain smaller 
sizes which do not drift quickly. Thus the high flux of 
material migrating into the inner regions from outside ex¬ 
ceeds the loss flux towards the sun, leading to significant 
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enhancements in residual solids. iCuzzi fc Zahnlel (|2004[1 
identified the same effect with respect to vapor, which 
can become enhanced inside EFs because its inward re¬ 
moval rate by gas advection is smaller than its supply 
rate, in the form of solids drifting from outer regions. 

It is illuminating to compare the mass enrich¬ 
ment/depletion results of the F, MTBF, and HA-MTBF 
cases (focusing on the red solid and dashed curves in the 
bottom panels of Figs. [Sj [7l and [151 the MTBF and 
BE cases are fairly similar). The most obvious differ¬ 
ence is in the outer nebula, where adding the bouncing 
barrier keeps particles, and associated radial drift loss, 
smaller. In case F, solids drain almost entirely from the 
outer nebula in 1 — 2 x 10® years, clearly problematic for 
outer planet formation. The HA-MTBF case represents 
another end of the spectrum, where the higher at keeps 
particles and their drift-loss smaller than the baseline F 
or MTBF cases, retaining outer nebula material longer 
and in higher abundance. 

In the innermost regions, both the F and HA-MTBF 
cases show the least buildup of material but for different 
reasons. In the F case, particles become sufficiently large 
that they can be removed by drift on short timescales, 
even while growing. This effect becomes very strong be¬ 
tween 10® and 2 x 10® years, after the outer nebula sup¬ 
ply has ended, and even the inner nebula starts to drain 
away. Contrast this with the HA-MTBF case in which 
particle growth to larger sizes is limited by strong tur¬ 
bulence, Stokes numbers are low, and thus so is particle 
drift. As we pointed out before, the temperature gradient 
is also very steep so that the gas advection velocity is also 
small. Thus the solids distribution changes much more 
slowly. In the outer regions, strong outflow transports 
mass outwards countering the inflow from drift of larger 
particles. These effects conspire to keep enhancement 
low, though a peak is beginning to grow at the water 
ice EF suggesting that at later times enhancements may 
increase. For both these cases, the solids enhancements 
never gets large eno ugh to lead to break through growth 
to larger sizes (as in iBrauer et al.l (|2008fl l. On the other 
hand, the slower particle growth rate of the MTBF case 
relative to the F case leads to smaller particles (but larger 
than the HA-MTBF case) which have slower radial drift 
rates, but relatively large enhancements of material in 
the inner nebula, at least at this stage in the evolution 
because of the strong mass inflow of material from the 
outer disk. 

Previous workers (e.g., iSteninski fc ValageasI 119961 
1997 : iQesla fc Cuzzil 200(I 1 g arau dl I2007I : iBrauer et al.l 


2008 : iBirnstiel et al.ll20ldr observed the transport of ma¬ 


terial from the outer disk to the inner regions, but gen¬ 
erally did not see the same total enhancements as we 
do here. In all these studies one or more simplifications 
were made, but it is probably still true that the differen¬ 
tial strength2»_^s_thejnosMmportantsm£le factor. In 
the case of iSteoinski fc ValageasI (|1996L I1997I1 . only icy 
material was treated, all solids vanished at the ice EF, 
no vapor tracking was included and a simplified particle 
growth algorithm was used, but some enhancement of 
solids in the ice-giant region was still seen due to inw ard 
drift from further out. While iBirnstiel et al.l (1201011 did 
study differential strengths, their models merely changed 
the particle strength at the iceline without removing the 
associated icy material, so it is unclear just how large 


an effect they would have predicted under more realistic 
assumptions. 

At the next level of detail, one sees discrete radial peaks 
and valleys in the abundance of solid material, some quite 
narrow and pronounced. These peaks generally coincide 
with locations just exterior to various evaporation fronts 
(EFs) due to water ice, and “organics”, which dominate 
the condensible material in our model. We suspect that 
for our fiducial case, if the silicate EF were completely 
resolved within our radial grid, that a peak may appear 
th ere as we l l. Sini ilar multiple features were pointed out 
bv iGaraiidI (|2007[ 1. The radial locations of these peaks 
vary with the radial location of the respective EFs, which 
are in turn dependent on the viscosity, temperature, and 
ultimately opacity of the disk. Thus, the system is cou¬ 
pled thermally and dynamically to the particle growth 
process. One can even discern fluctuations in the gas sur¬ 
face density associated with these structures (see figures 
Uli andUl and section WM- Enhanced local mass den¬ 
sity and/or smaller particles lead to higher opacity (as 
long as their size exceeds the peak thermal wavelengths), 
which affects midplane temperature and the EF. In case 
F, loss of solids leads to strongly falling opacity every¬ 
where, and a cooling nebula, with more rapidly evolving 
EF locations. In cases MTBF and HA-MTBF, EFs actu¬ 
ally evolve very little over the period of these simulations, 
allowing greater solids buildup behind them. The radial 
structure of the post-EF solids enhancements varies con¬ 
siderably, from an extremely narrow spike outside the wa¬ 
ter ice EF in case F (and MTBF), to somewhat broader 
“organic” buildups in cases MTBF and HA-MTBF (here 
a weak silicate buildup is present), each with associated 
opacity peaks. In fact, a glance at the red curves for Md 
in the lower panels of Figs. |4j|8l and fTHl shows that even 
at the last time depicted, convergen ce of solids mas s into 
narrow zones is still ongoing. While iGaraudI (|2007ll sug¬ 
gested that the amplitude of these peaks reaches a steady 
state in the long run, longer runs will be very instructive. 

4.2.1. Implications of mass redistribution and particle size 
for planetesimal formation: 

A subject of long-running interest is the still poorly- 
understood “primary accretion” of asteroids (in the 
inner nebula) and TNOs/KBOs/comet nuclei (in the 
outer nebula). Recently, attention has focused pri¬ 
mary accretion scenarios in which such objects “form 
big”, leapfrogging the various barriers to accretion 
of which only the first few are discussed here (for 
the others, see llda et al.l 120081: iNelson fc Gressell l201d : 

iGressel et al.ll2012l : lOrmel &: Okuzumill20I3ll . Two some¬ 

what disparate avenues to primary accretion may be 
distinguished: “streaming instabilities” (SI) which op- 
erate nearly axisYmmetrically in dense midplane layer s 
(iGoodman fc Pindc^ 1200(11 lYoudin fc Goodni^ 1200511 . 
and “turbulent concentration” (TC) which is a clumpier 
process, acting throughout the nebula but also lead¬ 
ing m ost effective l y to planetesimal s near the mid¬ 
plane (|Guzzi et al.l I2001L l2008l I20I0H . A full discus- 
sion of this top i c is b eyond the scope of this paper (see 
l.lohansen et al.l (j20I5ll for a recent reyiew). The condi¬ 
tions needed for SI are easier to assess, and we defer a 
discussion of the implications for TC to a future paper. 

For SI to operate requires a midplane particle-to-gas 
ratio Pd/P ^ 1 so the solids can driye the gas motions. 
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This has traditionally led SI advocates to focus on some 
combination of large particles, low turbulent intensity, 
and global enhancement of solids, usually due to late- 
stage removal of the gas. Related suggestions involve 
a reliance on “lucky par ticles” growing large enough to 
settle and lead to SI le.g.. l.Iohansen et al1l20I4fl, inaYbe if 
trapped at “pressure bumps” (iDrazkowska et al1l2013ll . 
but it seems to us that under nominal conditions, “lucky” 
particles (those of radius in our figures) are too rare 
and contain insufficient mass to drive collective effects 
such as SI (they might have other interesting implications 
as migrators, however). We believe it is more fruitful to 
explore the implications of the mass-dominant particles, 
represented by tm and StM- It is possible to express the 
local midplane value of ph/p, simplifying Eg. (fT4l a lso 
see Appendix B) as in iCuzzi fc Weidenschillina (1200(111 to 
hd/H = Y^oft/Stivi: 

- = ^^ = /|-=/\/si^. ( 62 ) 

p Zhd zj hd 

Figure 1201 shows Pd/P with other quantities as a func¬ 
tion of R (red dashed curve) for three representative 
models. In general, at this point in nebula evolution, 
Pd/P < (<C)0.1, but there is one intriguing situation 
in which pd/p may eventually approach a value where 
SI might be possible, and it is right where longstand¬ 
ing expectations are that the first giant planet would 
form - just outside the iceline (most prominently, up¬ 
per panel for LI). Moreover, this situation occurs quite 
early in nebula evolution and perhaps a longer simula¬ 
tion may eventually achieve this condition without the 
need to await photoevaporation or some other process 
to deplete the nebula gas. At this time and place, our 
results are showing particle radii of just over I cm, but 
recall that icy aggregates may not be able to compact 
fully and the identical StM could arise for porous par¬ 
ticles perhaps 10 or more times larger. Soon we will 
run our code for porous particles, and for longer times, 
to explore the possibilities for early planetesimal forma¬ 
tion. While SI appears inoperative in the inner nebula 
at these timescales for at values studied herein or larger, 
because the weak silicate aggregates do not grow large 
enough to settle into the required dense layers, the par¬ 
ticle sizes and density enhancements at these times are 
in the range of interest for TC, and also agree with what 
is seen in meteorites. It is increasingly suspected that 
some planetesimals did form in the inner nebula even at 
these early times, and it will naturally be of great inter¬ 
est to continue these simulations to the Myr timescales 
on which most observed chondrites form. 


4.2.2. Other implications of mass redistribution for 
observations of different kinds: 

The mass redistributions we see, fueled by the trans¬ 
port of material over large radial distances and the 
effects of Evaporation Fronts, have implications for a 
number of other currently puzzling observations. From 
the standpoint of mm-submm observations by ALMA 
and other radio observatories, the peaks and valleys in 
mass and opacity might provide a protoplanetary disk 
with a ringed appearance such as see n recently in HL 
Tau ([Partnership. ALMA et akl I20I5|1 and MWC 480 
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Figure 20. Radial profile of pd/p near the midplane is shown in 
the red dashed line. Streaming instabilities require pd/p ^ 1. It 
appears that, for some cases (specifically the upper panel) SI may 
at some point be achievable over a narrow range of radii, right 
outside the ice EF - a location where the first giant planet has long 
been thought to arise. Of course, the absolute location where this 
circumstance occurs will depend on model parameters. 

(|Oberg et al. 11201^ . Our simulations to date have not 
included EFs for the “supervolatiles” thought responsi¬ 
ble for the structure in MWC 480, but it would not be 
surprising to see a zoned or banded structure similar to 
that created by the ices, organics, and silicates we are 
currently modeling. 

From the standpoint of meteoritics, a number of im¬ 
plications are evident. The enhancement of vapor in the 
inner nebula, a combination of evaporated water ice and 
“organics”, is especially large in the LI case (factor of 
^ 8). Dependin g on the state of the evaporated “or¬ 
ganics” (see Sec. I4.5|l . whether CO, CH compounds, or 
higher hydrocarbons, the redox state of the inner nebula 
cannot help being influenced by this. 

Another current meteoritics puzzle involves the anoma¬ 
lously high ratios seen in most chondritic ma- 
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terials relative to the sun (|McKeegan et al.]l2008D , an en¬ 
richment that seems to appear some 1 —2 x 10® years after 
the e arliest “CAI” minerals condense, and to grow with 
time (lQavtonll200l : iSimon et alJl2012 ; iKita fc Ushikubol 

1201211 . Current thinking is that the heavy isotopes are 
liberated by a UV self-sh ielding process in the u pper lay¬ 
ers of the outer nebula (jLvons fc Yound [200511 . become 
bound in icy part icles, and drift inward s to evaporate in 
the inner nebula (|Yiirimoto et al.l[200lTH . Photochemical 
models seem to be able to produ ce the needed material in 
a timely way (iLvons et al.ll2009ll but there has never been 
a demonstration that it could be transported to the inner 
nebula on such short timescales. Now, it seems at least 
several plausible models indicate that substantial outer 
solar system material can indeed reach the inner solar 
system within 1 — 2 x 10® years. More detailed models, 
coupling these transport models with formation 

rates in the outer solar system, are in order. 

The peaks and valleys we see in the radial abundance 
of solids, even quite early on, may provide insight into 
the longstanding problem of the “ small Mars” relative to 
Earth and Venus (|Chamberil2001h : planetesimals formed 
from such material through the entire terrestrial planet 
region might have been quite nonuniform in initial sur¬ 
face mass density. Naturally, the absolute locations of 
peaks and valleys depend on the actual radial locations 
of the various EEs, as they evolve with time. It will be 
of interest to see what longer runs reveal, but the mere 
occurrence of mass peaks and valleys with up to order- 
of-magnitude density contrast in the terrestrial planet 
region on radial scales of an AU or so are surely of inter¬ 
est. Moreover, all of our models show large-scale inward 
drift of outer nebula solids, to varying degrees. This ef¬ 
fect alone might expl ain the compa ct surface solids mass 
densities inferred by iDeschl (|2007fl from early, compact 
giant planet configurations, even while the enveloping 
gas disk, with a certain amount of dust, might extend to 
much larger sizes. It is also well-known that “normal” 
(i.e. Kepler) planetary systems are much more packed 
with inner planets than is our solar system. The possi¬ 
bility of finding systematic causal factors for these im¬ 
portant differences is intriguing. 

4.3. Evaporation fronts; local refineries 

Just outside EEs, we find significant peaks in the solids 
mass, created when drifting, vaporized material diffuses 
back outside the EEs, recondenses on small grains, and 
advects or diffuses some distance before becoming accu¬ 
mulated into a larger particle. These peaks and valleys 
in the fractional mass of solids, which move as the EEs 
evolve, may be observable in protoplanetary disks, as 
mentioned in Sec. 14.21 Additionally, we see that EEs 
lead to a distillation effect whereby the eomposition of 
material residing some distance outside each EF is en¬ 
riched in the volatile that has just passed through that 
EF, and returned to recondense, leaving behind the more 
refractory material it had evaporated fro m, which con- 
tinues to grow and drift radially inwards (|Garaudl 120071 
previously referred to EEs as barriers to inward drift for 
the associated species). The red and magenta curves in 
the middle panels of Figs. El 0 and ITS] track this compo¬ 
sitional distillation through the internal density of local 
solids. For example, just outside the water ice EF the 
particle density drops from its “cosmic abundance” value 


of roughly 1.5 to a nearly pure water ice value of 1, indi¬ 
cating that the solids there have become highly enriched 
in water ice relative to cosmic abundance. Of course, the 
same effect would occur outside each EF and would occur 
for the various EFs we are not tracking here as well. It 
is a way of creating radial belts of solids with very differ¬ 
ent composition than material on either side. The radial 
width and shape of these refinery or distillation zones de¬ 
pends on at primarily, but to some degree on the nature 
of the growth-by-sticking assumptions. Moreover, the ra¬ 
dial zones evolve with time as the EFs and the mass flux 
through them vary, so longer simulations will be needed 
to what happens at later times. One can imagine inter¬ 
esting meteoritic and astrophysical chemistry associated 
with the EFs for silicates, troilite, and organics as well. 
This is a potentially rich area which also has implica¬ 
tions for circumplanetary disks, and will be developed in 
future studies. 

One prediction can be made in this regard, which we 
should mention here. The ice-enriched nature of solids 
immediately outside the ice EF, combined with our pre¬ 
ceding results mentioned above, that the most likely 
place for early planetesimal formation to occur by SI 
is just outside the water EF, leads us to predict that 
Jupiter, say, if that were the planet first formed just out¬ 
side the snow line, wherever that arose due to energetics 
details, would have a water-enhanced core relative to cos¬ 
mic abundance. Perhaps Juno can test this prediction. 

4.4. Nebula opacity 

We have followed the evolution of EFs using a self- 
consistent calculation for the photospheric and mid¬ 
plane temperatures, that employs a new model to com- 
pute the Rosseland mean opacity of aggregate particles 
(|Cuzzi et al.l 12014'^ . The opacity depends on particle 
size, number density and composition, so it depends on 
the evolving particle size distribution and abundance. 
We have found that the opacity in the outer disk be¬ 
comes quite small if growth proceeds quickly (particles 
growing larger than a thermal wavelength have smaller 
cross-section per unit mass), and/or if the loss of mate¬ 
rial to inward drift is high (less material). This tendency 
is weaker for models with higher at that keep particles 
smaller, allowing them to be retained longer and be more 
easily carried outwards as the disk spreads. In the inner 
region, the opacity can become larger with time, even 
though particles are growing to larger sizes (although 
slowly) if supply of material from the outer regions can 
keep the mass densities high. A sustained high opacity 
can maintain higher temperatures in the inner regions 
for longer periods. We also note that the radially vari¬ 
able opacity can even, apparently, lead to “lumps” in the 
gas surface density, because the opacity affects the local 
temperature and thus the viscosity in our a-model. 

The evolution of opacity is important for two new hy¬ 
drodynamic instabilities which apparently lead to sus¬ 
tained turbulence, without any need for magnetic fields 
or ionization, and are thus extremely important to the 
inner solar system which has been regarded as “dead” 
or at best “ d ull” fr om the stan d point of turbulence 
(iNelson et al.l [201^ [Turner et al.l 12014 iStoll fc Kiev I 
12014 IMarcus et al.M2013l 120141 . The new instabilities 
(as currently understood) require a suitably rapid ther¬ 
mal equilibration, and thus prefer lower Rosseland mean 
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op acities. The l imitin g opacities have been estimated 
bv iNelson et al.l (1201311 to be on the order of 1 cm^ 
(perhaps up to 5 cm^ g“^; Umurhan, personal communi¬ 
cation 2015). In the outer nebula, our opacities remain 
at or below these nominal thresholds, but in the mass- 
enhanced inner nebula regions, opacities reach larger val¬ 
ues. Whether the corresponding instabilities are damped 
or not under these conditions is thus a pressing topic for 
future study, perhaps even in a coupled fashion with par¬ 
ticle growth. 

4.5. Caveats and Future Work 

There are several areas in which the code can be im¬ 
proved, which will be introduced in later papers. The 
initial conditions for the models in this paper are consis¬ 
tent with the earliest times in the nebula, but we have not 
added ongoing infall from the parent cloud. Such effects 
can b e included through relatively simplified ap proaches 
(e.g.. lBirnstiel et al.ll2010HYang fc Cieslall2012ll . 

None of the models presented here have reached a sit¬ 
uation in which the mass fraction of vapor in the en¬ 
hanced inner nebula exceeds the local H 2 -He gas density, 
but at least one has come within a factor of ^ 10 (the LI 
model. Sec. 13.6p prompting us to consider that this situ¬ 
ation may arise for a given model or at longer simulation 
times. Since this added vapor affects the disk viscosity, 
it is important to include its effects when this occurs. 
Furthermore, the local gas density influences the process 
of grain growth through the relative and radial veloci¬ 
ties, thus a more self-consistent treatment requires that 
the volatilized-condensible vapor component is included 
in any calculations that involve E. 

There are several areas in which we can improve the 
treatment of particle growth and composition. Most 
pressing is to explore the role of porosity. The moments 
method we utilize for co agulation is ideally suited for this 
(|Estrada fc Cuzzill2008ll . The fractal growth of particles 
by low velocity sticking of small monomers will cause 
the particles to have a much lower densit y than the solid 
densi t y of the monomers t hemselves (e . g..lDominik et all 
120071 : lOrmel et al.1 120081 : iZsom et 1201011 . Coupled 
with the stickiness and strength of icy particles, it has 
been argued by some workers that these fluffy aggre¬ 
gates can grow large, but maintain low relative veloci- 
ties, allowing for cont i nued growth without com paction 
(jOkuzumi et. al.1 120121 : lOrmel fc Okuzuniil I2013I1 . In the 
extreme case that a particle grows with a fractal dimen¬ 
sion of exactly two, the stopping times remain the same 
for a single monomer. Thus, a model like our Q107 case 
with higher porosity might allow particles to grow faster 
and overcome the radial drift barrier. Furthermore, as 
we have noted in Sec. 12.4.41 our assumption of a gaussian 
PDF for the relative velocities in our calculation of the 
migrator destruction probabilities has b een called into 
question from more recent w ork (e.g., iGalvagni et al.l 
1201 It iWindmark et al.ll 2012 bll. thus we intend to explore 
the more rigorous model of iGaraud et al.l (1201311 in fu¬ 
ture modeling efforts. The choice of PDF may very well 
affect the frequency of lucky particles which may modify 
some of our conclusions here. 

Other improvements concern our treatment of EFs. 
Currently we make the simple assumption that inwardly 
migrating particles lose all their volatiles at the corre¬ 
sponding EFs, while modeling the EF as a radial range 


of nearly constant temperature, centered on the species 
evaporation temperature and buffered by evaporation of 
the associated source of opacity. Though this works well 
for particles of the size we are so far growing, it may 
fail to capture the case of a much larger particle where a 
particle’s surface layers can insulate its material at depth 
while it drifts a considerable distance. Further sweepup 
of local refractories may effectively trap any volatiles re¬ 
maining deep within the particle and offer a way to de¬ 
liver volatiles deeper into th e inner regions of the nebula 
(e.g., iSuDulver fc Linl[200(Tll . 

We have also assumed that evaporation and conden¬ 
sation are reversible processes for simplicity, but this is 
certainly not always the case. In particular, the “evap¬ 
oration” of “organics” may instead be a decomposition 
of the original complex molecules into nearly inconden¬ 
sible species like CO or CH 4 , which participate in gas- 
phase chemistry rather than recondensing outside of the 
“organics” EF. This loss of refractory, carbon-rich or¬ 
ganics must, at some level, be a one-way street, because 
none of the most primitive carbonaceous chondrites con¬ 
tains more than about one-tenth of the carbon found in 
comets (see Cuzzi et al 2003 for a discussion and refer¬ 
ences). This loss of refractory organics may effect par¬ 
ticle strengths as well, because organics are often found 
to be “sticky”, perhaps extending the range of “strong” 
particles to locations somewhat further inwards than the 
water ice EF. 

Our radiative transfer treatment (section [2)3]) assumes 
implicitly that the bulk of the energy released by vis¬ 
cous angular momentum transport is dissipated near the 
midplane, or at least in a fashion proportional to gas den¬ 
sity, and is transferred vertically upward and outward to 
space, making the midplane the highest temperature lo- 
cation. However, in cases where visc osity is MRI-related 
(iTurner et al.ll201" il. or VSI-related (jNelson et al.ll201^ 
iStoll fc Kiev 112014 ). it may be that the dissipation max¬ 
imizes at some altitude near one scale height, or even 
closer to the disk photosphere. Such a case can surely 
be modeled and probably parametrized for our evolu¬ 
tions. The magnitude of the effect can be estimated by 
considering that, in the worst case where the energy is 
dissipated right at the photosphere, half of it is lost to 
space immediately and half radiates downwards. So, the 
vertical temperature distribution (that we do not model) 
will become more nearly isothermal, but the midplane 
temperature at a given location may not change by more 
than a factor of 2 ^/'*. 

We expect there will be other innovations as our un¬ 
derstanding of the problem improves, but our code is at 
least well suited and robust enough to readily accept all 
the new physics mentioned above. An additional layer 
of parallelization in mass bins will greatly speed up our 
code which will allow for more complexity, and will be 
implemented in the future. Finally, we note that this 
model can readily be applied to the circumplanetary disk 
environment. 
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APPENDIX A: CODE ALGORITHMS 

In Sections A.1-A.2, we describe how we separately 
solve the equations for the radial evolution of the gas 
(Eq. A-[Tj), and that of the vapor and dust components 
(Eq. A-[TT]) in the first two sections. In Sec. A.3., we 
then describe in more detail how the temperature and 
opacity are calculated. The calculations in Sections A.l- 
A.3 are all carried out at synchronization steps and thus 
are only done at regular intervals. In Sec. A.4. we fur¬ 
ther describe the details of how we solve for the relative 
velocities which are required for both coagulative and 
incremental growth. Sections A.5 and A.6 describe the 
algorithm for growth beyond the fragmentation barrier, 
and radial drift of migrators. 

A.l. Gas Evolution Equation 

We use an implicit finite differencing scheme to solve 
Eq. (III). Radial bins j are spaced logarithmically over 
the radial range Rin < Rj < Rout such that 

Rj = Rin(Rout/Rin) , (A-1) 

with Ri = Rout and Rn^ = Rin, and riR, is the total 
number of radial bins. We thus first rewrite Eq. © 
in logarithmic form, using equation and the relation 
i' = OtcR = Otc^/ri: 


'm 


3A d 

^51nR 




(A-2) 


where A = Finite differencing of Eq. 

(|A-2|) leads to the tridiagonal system -|- 6jEJ^^ -|- 

where n is the time 
index for time t and n -|- 1 for time t + Atsync- The 
coefficients are given by 
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(A-4) 

where Ej’ = ctJ’^E, the first term on the RHS side is the 
diffusion term with ZIv,d the diffusion coefficient, and the 
second term on the RHS is the advective term. Dropping 
the specifics of whether we are treating vapor or the dust 
for now, and also suppressing the species subscript i, 
finite differencing again leads to a tridiagonal system for 
each species i, where the coefficients are for the diffusion 
term (1st term on RHS) 


ia‘E,r,R-i/ 2 AlnR(Ad^i - A^, , 

ia‘E,R,R-i/ 2 AlnR(Ad^i - A^,,)^, 

(A-5) 

with S = (l/2)Atsync/(A InR)^, and now, recalling Eq. 
(dH) for Dd, 




Pkj A 
1 -I- Sti j 


(A-6) 


for the dust component, whereas A^ = A for the vapor. 
The advection term (2nd term on RHS) requires that 
we use upwind/downwind differencing for the derivative 
that depends on the sign of the velocity 14,d- If 14.d < 0, 
we add terms to the Uj and bj: 


a, = a, - <5AlnRI^4VAl,-iE,_i/R|_i/2, ,, 

bj = bj + 5AlnRI/”+iRjEj/R2_^/2, ^ 

whereas if 14,d > 0, we add terms to the bj and Cj'. 


Qj — 6AR^_y^{Rj_i/Rj)‘^a^j_^Tj-i, 

b, = l + SAiR-y^%+R-li%)a)T„ (A-3) 

Cj = —5AR^.^(^y2(Rt+i/Rj)^Q^5'+iTj+i, 

with Cj = —ttj, gj = —Cj, fj = 2 — bj, and <5 = 
(3/2)A4ync/(AlnR)^. Since bin spacing in our code is 
logarithmic, the j ± 1/2 refer to the logarithmic center 
between bins j and j ± 1. 

A.2. Advection-Diffusion Equation 

The advection-diffusion equation can be differenced in 
a similar way to that of the gas evolution. We begin by 
rewriting the RHS of Eq. m in logarithmic form 


— ^3 dAhiRVj ’’’ RjT,j/R-j^-^!^, (ts,Q\ 

c, = c, + 5AInRI/"+V^.+ iS,+i/R 2 +i/ 2 - ^ ^ ^ 

Finally, the remaining coefficients are Cj = —Oj, gj = —Cj 
and fj=2 — bj with velocities evaluated at current step 
n. 

For the inner and outer edges of the disk, we impose a 
combination of Dirichlet and Neumann boundary condi¬ 
tions (Robin boundary condition) through the mass flux 
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where is the mass flux of vapor (or solids) of species 

i. The boundary conditions are also made implicit, and 
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require the presence of ghost points just outside the 
boundaries. We extrapolate to find the values of the 
necessary quantities using polynomial methods. 

Application of the inner (j = tir) and outer {j = 1) 
boundary conditions requires that w e alt er the tridiago¬ 
nal coefficients. We difference Eq. (IA-91) by taking the 
derivative centered around the boundary points. In order 
to estimate the flux at time for species i, we define 
/fj’’^ where is the total mass of 
solids or vapor in a bin. The modified coefficients for 
the boundaries, again dropping specifics of dust or vapor 
with the exception of the total fractional masses fj’’^ to 
avoid confusion, are then 


b\ = 


bi-ai ( 2 Vj 


n +1 


fi = fi-ei 


AlnR 




Mi/7ri?i/rS? ; 
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AlnR 
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(A-IO) 


Aggregate particles and Effective Medium Theory 
(EMT). We assume in our model that particles are com- 
positionally heterogeneous, being granular aggregates of 
much smaller grains of all compositions that condense 
at the local temperature, as is consistent with obser¬ 
vations of meteorites and interplanetary and cometary 
dust particles. Assuming our particles to be aggregates 
composed of much smaller elements which are them¬ 
selves smaller than any relevant wavelength allows us to 
use Garnett Effective Medium Theory to calculate av- 
era ge refractive indice s of their ensemble (see Appendix 
C, Cuzzi et al.l l2014all. The Garnett model (see, e.g., 
iBohren fc Huffmanll 19831) describes inclusions of some di¬ 
electric constant cq embedded in a homogeneous medium 
of dielectric constant Cm (which is vacuum here), and 
can be generalized to a multiple component, potentially 
porous medium composed of multiple species with differ¬ 
ent dielectric constants and volume fractions relative to 
the aggregate volume. 

The monochromatic opacity kx is then calculated from 
our averaged indices (obtained using Garnett theory) 
with 
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KX = - 
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7rr^nd(T) [Qe(r, A)] dr k. 
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with c[ = ai + Cl, g'l = ei + gi, = a„R -I- c„r and 

~ SnR + ffnR- 


where we have approximated the integral by a summa¬ 
tion over all particle sizes k weighted by their discretized 
number densities bk (cm“^), and the are Simpson’s 
trapezoidal rule coefficients. Note that nd{r) dr has been 
replaced with the bulk number density bk which applies 
to both the dust and migrator populations. The bulk 
number density for the dust population is given by 


A. 3. Opacity 

We employ the opacity model of IGuzzi et al.l (|2014aD 
in order to calculate the variable Rosseland mean opacity 
K of the evolving particle size distribution. Galculation 
of K involves determining the wavelength-dependent, or 
monochromatic opacities kx for a gas-particle mixture, 
which depend on an extinction efficiency Qe = Qs. + 
Qs for grain size r and wavelength A. The extinction 
is due to a combination of pure absorption, which is re¬ 
radiated as heat, at an efficiency Q„ , and sc attering at an 
efficiency Qs- Following iCuzzi et al.l (l2014all . we allow for 
nonis otropic scattering b y adopting the common scaling 
(e.g., Ivan de Hulstlll980l) 

Q's = Qe-gQs = Qn + il-g)Qs, (A-ll) 

where the asymmetry parameter g = (cos 0 ) is the first 
moment of the distribution of the scattered component, 
or phase function P{&), w hich describe s the de gree of for¬ 
ward scattering {g > 0, see lGuzzi et al.l (I2014a() . their Eq. 
3). Values of 5 < 0 are preferentially ba ck-scattering and 
g = 0 is isotropic. IGuzzi et all (|2014aD give expressions 
for the absorption and sc attering efficiencies in terms of 
the electric polarizability (|van de HulsdUfiSOl) . which is a 
function of the particle size and complex particle dielec¬ 
tric constant. 
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In Eq. (IA-I3I) . pd is the volume mass density of total 
solids in the column, mmin and wl the minimum and 
maximum masses in the dust particle size distribution 
(Sec. 12.4.IL and subsequently bk = b’f for migrator par¬ 
ticles (Sec. 12.4.41) . The Rosseland mean opacity is then 
calculated from Eq. (ED). 

The algorithm for calculation of the midplane temper¬ 
ature for each radial bin Rj is called at every tsync in 
conjunction with the gas evolution and solids and va¬ 
por advection and diffusion. The calculation is iterative 
because the layer optical depth depends on the opacity 
through the evolving size distribution, which depends on 
Tj. We use a root finding technique to evaluate Eq. ED 
in which each temperature estimate requires that we re¬ 
calculate the dust fractions af j and particle densities p^, 
taking into account that the radial location can be an EF 
(Sec. 12.4.51) . to determine k. The new opacity estimate 
then determines a new Tj. Iterations continue until the 
change in Tj on any subsequent iteration is less than 
some tolerance level. 

During the iteration process, there are as many as five 
possible circumstances regarding the previous bin tem¬ 
perature Tj and the new temperature Tj, 
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1) Tj < T, + ATef, T' > T* + ATef; 

2) Tj >T, + ATef, T' <T, + ATef; 

3) Tj, T' <Ti + ATef; (A--14) 

4) Tj > Ti+i — ATef, Tj < T^+i — ATef; 

5) T, + ATef < T,,Tj < T' < T,+i - ATef; 


which we consider when adjusting a radial bin’s solids 
and vapor content. For example, in condition (1) and 
also condition (3) for the case in which Tj < Tj^ changes 
involve the conversion of solids of constituent i into va¬ 
por: 


a, 


df 

m/ 

V/ 


= af^ (T, + ATer - Tj)(T, + ATer - T,)-i; 

= a-(T, + ATer - Tj)(T, + ATer - T,)-i; 

= + ^<4,0 + 

(A-15) 


where, e.g., ^4,j — 4,j ~ 4[j- the case in which 
all remaining solids are converted to vapor (condition 
1) then afj and a™ are simply added to a^ j. On the 
other hand, in a situation where vapor is condensing, for 
example condition (2), or also (3) but with Tj > Tj then 
we partition vapor onto both dust and migrators in the 
appropriate area fractions 


= al^{T'^ - T, + ATer)(T, - T, + ATer)"!; 

afj = 4 j + x^jAaJj-, 

< = + (1 - 4)^4.v 

where Xs is the area fraction of dust in solids. The it¬ 
erated values of the mass fractions are then used to cal¬ 
culate the mean particle material densities and adjust 
the particle mass distributions. The final values of the 
mass fractions are assigned when the temperature has 
converged. 

A special case in our code involves the conversion of 
troilite. Solid troilite is converted to vapor and metal¬ 
lic iron in the fraction fi = 56/88. For simplicity, we 
consider this reaction to be completely reversible, where 
more care must be taken in the reverse process to assure 
the fractions are correct. 

A.4- Relative Veloeities 

In our calculations, we assign a single particle-to-gas 
and particle-to-particle relative velocity for the entire 
particle mass distribution (dust and migrators) within 
the subdisk “dust” layer, which is defined by the scale 
height /iD as described in section 12.2.31 This calcula¬ 
tion uses a logarithmic spacing in particle radius from 
the minimum size rmin which is fixed in our code, to the 
largest size ru. 

The stopping times of solid particles t% depend on the 
flow regime which depends on whether the radius of a 
particle is larger or smaller than the gas molecular mean 
free path (Sec. I2.2.1|) . For any height z, the stopping 
time for particle of mass rrik at radial location Rj is 


4 = ^4pt/ T^rlpt + E 
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where the gas density as a function of height is pj = 
Y^2/7r(Ej/2iJj ) exp [—i(z/i7j)^] and the expression in 
brack ets is the Epstein-to- Stokes regime transitional for¬ 
mula (|Podolak et al.l[r988f) . with coefficients At = 1.249, 


Bt = 0.42, Ct = 0.87 and 


Dt,j = -(1/e) [At + Bte-^4 . (A-18) 

where t'm = Pm/P is the molecular viscosity, and = 
1.3 X 10“"'^ g cm“^ s“^. The transition between the 
Epstein and Stokes regimes is defined by < eA™^'’ 
where the factor e = 9/4 is typically used. In our code 
we have adjusted this value to get a smoother transi¬ 
tion and found 3/2 works better. For Vk A™^^, the 
stopping time formula reduces to the Epstein regime 
tl = TkP^CjPj- For Tk > utilize the Stokes 

regime equation 
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where AV)f® is the particle-to-gas relative velocity. The 
drag coefficient Cj) depends on the particle Reynolds 
number RePfc = 2rfcAl^P^/^'m and i s calculated using the 
prescription (|Weidenschilhng| [197711 : 


p4/RePfe forRePfc<l; 

= < 24/ReP° ® for RePfc < 800; . (A-20) 

[0-44 forRePfc> 800; 

Since depends on AV)f® for RePfc > 1, iterations will 
be required to determine the correct stopping time, and 
thus relative velocities. 

We begin with equations (I34II37I1 and solve for the gas 
radial and azimuthal velocities u and v at semimajor axis 
Rj as 


Efc Akpt [BUk + 21114 ] + 2BnpV^ 

(A-21) 

B^ + U? 

- Efc Akpt [{n/2)Uk - BVk] - n^VK 

(A-22) 

52 + 


where B — Aipf. These can then be inserted into the 
equations for the Uk and 14 and the individual compo¬ 
nent can be removed from the sums. That is, we solve for 
the component 14 (or f7fc) by substituting the expression 
for V into the component equation for 14 (or Uk, e.g., Eq. 
[33] or [33]). This gives 


{n/2)FkUk + Yk+n^Akpr,VK 
AkpGk 


(A-23) 
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where Gk = BAkp'l - B'^ - Fk = A^pp'^. + 
and 


n = AkP Y, Aipf - BVi]. (A-24) 

l^k 

Solving for the final expressions for the velocity compo¬ 
nents, we find 


Uk 


AkpGkZk + 2QFkYk + 2QAkp{i^^Fk + 


AkpBGk)pVK 


_ , (A-25) 

[{AkpGk)^ + {nFkf] ' 


Vk = - 


{p,/2)FkZk — AkpGkYk + Akpfl^{BFk — 


AkpGk)pY^K 
where Zk is given by 


[{AkpGkf + {nFkf 


(A-26) 


Zk = Akp Y \BUi + 2m] ■ (A-27) 

l^k 

The equations for Uk and Vk represent 2n equations that 
can be solved by matrix inversion using standard meth¬ 
ods. 


A.5. Migrator Growth and Radial Drift 

If we define the mass of the dust subdisk in a radial bin 
Rj to be = 2pjhj£Yj where the area of a radial bin 
is jYj = 7r(i?^ — and the subdisk layer height is 

defined in Sec. I2.2.3p . then the total mass of new “mi¬ 
grator” material crea ted in time Atj is roughly (compare 
to Eq. of section 1^.4.5|) : 


M* ~ Mf 




2-9 \ 

— m* ^ \ 

2^ I > 

- "*min / j 


(A-28) 


Although the migrator population covers masses m > 
niif, we consider the mass in this layer including 

the fraction of dust lying in this layer, to be available 
for sweepup by all migrators. This mass does not reflect 
all the mass in radial bin j, as some fraction of “dust” 
occupies the upper layers beyond /id- The migrator mass 
thus “created” in time Atj is assigned to the first total 
mass bin of the migrator population, = M*, while 
already existing total mass bins are shifted up by one in 
index. 

Drift mass. Because we use an asynchronous time step 
scheme in our code, the timestep of a radial bin j is typi¬ 
cally longer than that of bin j — 1 which lies exterior to j. 
In order to account for the asynchronicity, we use a pre¬ 
dictor/corrector method for the migrator mass drifting 
in. There are three circumstances that arise within this 
scheme: (1) a time step for bin j — 1 was called at the 
same time as bin j; (2) no time step was called for bin 
j — 1 when a step for j was called; and, (3) it is a synchro¬ 
nization step in which all bins are called, but with the 


appropriate time AtY'^'^ set for each radial bin such that 
all radial bins are brought to the same elapsed time. The 
synchronization step is the step in which all drift mass is 
accounted for. The three cases can be summarized as 


Mpc = 


m: 


drift 

i-1 


- {Nj-1 - I)Atj_i 

At,., 

,N-i [Nj - l)Atj - {Nj.i - l)Atj-i 


Mff:f^{At,/At,.,y, 


Atj. I 


^drift,A _ ^drift,A-l^ 


{Nj - l)AtY - {Nj.i - l)AtYf 
A^ ’ 

(A-29) 

where = NgyncAt^in — NjAtj, and the superscript 

N refers to the most recent stored calculation, and N —1 
the calculation previous to N. 

Due to the nature of our parallelization of the code 
in radial bins, the most recent information for the drift 
mass is not available until a synchronization step. An 
individual radial bin j is assigned to a single processor, 
but requires information from j — 1 regarding the mass 
drifting in. We store and use values from previous time 
steps for the mass that drifts out of bin j—1 into j 
and we use those as the basis of our calculation which 
leads to some small “out of phase” transport which is 
eventually corrected at a tsync- Given that the amount 
of mass drifting in compared to the total mass of a bin is 
quite small, any temporary effect that may occur during 
regular time steps is minimal. 

Afpc represents a total mass of migrators of various 
sizes available to drift inwards into radial bin j in time 
Atj.i, whose total mass per mass bin irik is given by 
such that The actual mass 

available in the time step Atj (or in the case of a syn¬ 
chronization step, AtY’^^) is then scaled to M^Y-i = 
{Mpc/Mj!ji)M^^Y-i- The particle mass distribution for 
bin j — 1 will be different from bin j, thus we reinterpolate 
the drift mass distribution from bin j — 1 onto the same 
particle mass range as bin j and add the contribution to 
the total migrator mass distribution M^Y ■ lu the event 
that the particle mass range in bin j — 1 does not all lie 
within the mass range of bin j, we (1) return to the dust 
population of radial bin j the mass of all migrators that 
are in sizes smaller than the fragmentation barrier mass 
of bin j; and/or, (2) add new particle mass and total 
migrator mass bins to radial bin j for particles from bin 
J — 1 that are larger than the most massive pre-existing 
particle in bin j. 

Evaporation fronts. If radial bin j or j — 1 (or both) lie 
within an EF, then inward drifting migrators will have 
their compositions altered through adjustment to the lo¬ 
cal nebula temperature. This means that some or all of 
their volatile will be removed and added to the vapor 
content. For an EF with evaporation temperature R, we 
adjust the mass fraction in species i of migrators from 
bin J — 1 
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(tef = 


f Ti + ATef — Tj \ 
+ ATef — Tj-i) 


(A-30) 


with the evaporated fraction added to the vapor content 
of bin j : 


drift/ 




^i-1 -^j-i ^j-1 


£/j 


M, 


(A-31) 


t-i 


In Eq. (IA-311) above, is the adjusted total mass 

of migrators drifting in to bin j in time Atj, is 

the total mass of migrators in bin j — 1, the ratio of 
which along with the ratios of gas density and bin area 
ensures we are adding the correct mass fraction to bin j. 
In the special case of the troilite EF, the above equation 
will have a (1 — //) factor multiplied to the second term 
on the RHS, and we add the remaining factor // to the 
metallic iron dust fraction (see Appendix A.3). 

Migrator growth. The incremental growth of migrators 
is described by Eq. (l4^ . We assume that all migrator 
particles rrik.j grow from both the dust content in the 
subdisk layer, and potentially other migrators to some 
new size m'j. ^ in time Atj. The fraction of dust and mi¬ 
grators accreted can be determined from the individual 
terms on the RHS of Eq. (1441) . The change in the total 
mass of migrators in mass bin k is estimated then from 


fraction over all masses accreted by I , and is thus related 
to the total mass of migrator material available in the 
subdisk for a migrator to accrete from. The factor fp is 
the fraction of dust accreted in Atj. 

If it is found that more mass from a migrator bin is 
accreted than what is available (be., M™' < 0), we define 
a scaling coefficient for each bin k found to have more 
mass accreted from it than available: 


-4/cj- = Ml 


k,j 


-1 


(A-34) 

which weights the amount of material that is subtracted 
from a bin due to accretion from other bins, and be¬ 
come the coefficient for the sum on the RHS of Eq. 
(IA-331) . The overestimated mass is then l — Akj which is 
subtracted out of AM^^ once iterations that ensure the 
modified condition for Eq. (jA-331) are satisfied for all k. 

Once all growth calculations are done, as a final step 
we perform a reinterpolation of the mass histogram of mi¬ 
grators (including all related quantities) to a logarithmic 
mass bin spacing, and make sure that the criterion for 
the maximum number of bins per mass or radius decade 
is adhered to (see Drazkowska et al. 2013, 2014). 


A.6. Particle Destruction Probability 


AM“ = (- 1 , (A-32) 

\mk,j J 

which can be an overestimate for two reasons. 

The first reason is that there is a limited amount of 
dust mass available in the subdisk that can be accreted. 
If it is found that more dust mass was accreted than lies 
within the subdisk, the excess is subtracted out and the 
particle sizes, densities and accreted mass are adjusted 
accordingly with the same factor taken from all sizes. Al¬ 
though the subdisk dust mass becomes negligible, more 
dust remains at higher altitudes, which can be mixed 
downwards by settling and diffusion. 

The second reason is that the initial calculation can 
allow more mass to be accreted from migrator mass bin 
k by other migrators than is available in mass bin k. 
Under these circumstances, the entire redistribution of 
the migrators will need to be redone (perhaps multiple 
times) until only the mass available for accretion can be 
accreted. Under normal circumstances, after the growth 
calculation, we redistribute total migrator mass accord¬ 
ing to the fraction of Eq. (IA-321) that was due to mi¬ 
grator accretion of other migrators. For every migrator 
total mass bin fc, we loop over other migrator total mass 
bins I > k (here I and I' represent dummy indices) and 
subtract out what those bins accreted from mass bin k 


C 

MZ = - (1 - /p) E 


where Sk.ip’kj represents the density fraction of mate¬ 
rial of particle mass bin k in the migrator population 
accreted by mass 1, and is tii® total density 


Rather than adding another time consuming calcula¬ 
tion in our migrator growth subroutine that would in¬ 
volve integration over multiple integrals to properly cal¬ 
culate the destruction probabilities, we discretize the de¬ 
struction rates of a migrator particle of mass nik due to 
dust and to other migrators as 


= (A-35) 

/=1 


k 

= E (a-36) 

i=i*+i 


where summation to Z = /» is over all dust masses up to 
the fragmentation barrier to*, and the summation from 
/* -|-1 to fc is over all migrators of mass rrik and smaller. 
Equation (IA-36P for migrators contains an additional fac¬ 
tor that depends on the probability of destruction of 
migrators during the current time step. We account for 
the second integral (over time of growth) by using a dis¬ 
crete number of points over the time interval as it grows 
from TOfc to to).. Different mass values between mk and 
to)., along with their associated relative velocities, deter¬ 
mine the destruction probability over the timestep. Be¬ 
cause growth is relatively small in time At^- over a wide 
range of values of turbulent intensity at , a few points is 
generally enough, but more points are used when growth 
is fast, such as it might be the case for low at. The total 
probability of destruction of a migrator in mass bin k is 
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k = 


Ai 






1 — exp 


M 


p=i \i=i 

N„ / 


-irEI E 


p—1 \Z=Z* + 1 


(A-37) 


where the index p and Np refer to the number of points 
used in the estimate. Once I^k is determined, we can cal¬ 
culate the total mass of particles of size k destroyed as 
^dest _ this mass is returned to the dust pop¬ 

ulation. Bins that suffer total destruction are removed if 
they are at the tail end of the distribution only. Bins with 
zero mass in the middle of the distribution are allowed. 


APPENDIX B: PARTICLE DIFFUSIVITY AND SCHMIDT 
NUMBER 

A quick survey of the fluid dynamics literature reveals 
that the term “Schmidt Number” (Sc) is widely used 
as the ratio of gas viscosity v to particle diffusivity 
for particles of arbitrary Minor complications arise if 
Scg = vlD^ , for the gas its e lf, is significantly different 
from unity (IStevenso M ll99Ct iHughes fc Armitag^ 120101 
I 2012 f) because many deriva tions of Da relate particle dif¬ 
fusiv i ty to gas diffusivity (iVolk et al.l 119801: iCuzzi et ahl 
1993t iCuzzi fc HogjanI 12601 ^ iSchraoler fc Henning! 12004 
Youdin fc Lithwick 1200711 we will assume Scg = 1 , but 
deviations are easily allowed for. He re we have adopted 
the fo rmulation of and Sc in lYoudin fc Lithwick! 
(|2007L henceforth YL07), who incorporated orbital dy¬ 
namics effects . For discussion and vali dation of these for¬ 
mulations see iCarballido et all (|2011ll . who also showed 
that the resulting particle layer thickness varied as 
St“^/^ for St = 0.1 to 100. This innoc ent-looking St- 
depen dence, familiar from the results of iDubrulle et all 
(119951 henceforth D95), hides several subtleties, and is 
not quite the final answer, as we sketch below. 

A closed-form expression for the particle vertical scale 
height hd can be obtained by integrating the verti¬ 
cal mass balance equation for the concentration ratio 
Pd/p (D95 Eq. 32; Garaud et al. 2007, henceforth 
G07, her Eq. 21 ), or treating mass balance us ing scal¬ 
ing arguments (IGuzzi fc Weidenschillingl I2006L hence¬ 
forth GW06 (their Eqns. 4-7)). Written in terms of 
Sc, without worrying about its proper form (see however 
ISchrapler fc Henning! (I2004f) for a pre-orbit-dynamics as¬ 
sessment), Eq. (22) of G07 for hd/H is exactly repro¬ 
duced (ignoring constants of order unity) by combining 
Eqns. (37) and (39) of D95 (but note that a factor of 
= 1/Sc has been lost in D95 Eq. (37), because Eq. 
(36) requires the normalized diffusivity Kt = Ktocq/St). 
Solving the diffusion equation for a particle/gas ratio is 
a preferred approach; because neither CW06 nor YL07 
accounted for vertical variation of gas density; their ex¬ 
pressions give hjd > H when St < at, which is unphysi¬ 
cal. 

GW06 assumed a constant gas density and so their 
result is directly comparable to the scale height h/H for 
the ratio Pd/p derived by D95 (her Eq. 36); in fact, it is 


identical if both are expressed in terms of Sc generically. 
Thus, it can also be transformed to a proper particle layer 
scale height hd/H using Eq. (39) of D95, also becoming 
identical (in terms of Sc, for a single particle) to Eq. (22) 
of G07 (again neglecting constants of order unity): 

hd/H = (l-kStSc/at)”^/^ (B-1) 

However, all these previous results from mass balance 
equations assumed a particle “settling time” fsett given 
by the terminal velocity of a single particle under ver¬ 
tical solar gravity = l/StH). This is not appropriate 
for St ^ 1 particles whi c h are under damped, as pointed 
out by iWeidenschillind (|1984 1198811 who suggested us¬ 
ing the inclination damping time of a layer of large 
particles (= ts) as the “settling time” when St;^ 1 . 
T his idea was incorp orated into the numerical models 
of IGuzzi et al.l (1199311 and used in a bridging expression 
constructed by YL07 (their Eq. 2) to derive their Eq. ( 6 ) 
for h/H, which ironically reduced exactly to the small- 
St (Sc = 1) expression in CW06 or Eq. (37) of D95. 
This is because (relative to the earlier expressions), the 
larger Sc = 1 -I- St^ of YL07 (using their physically mo¬ 
tivated Eq. 5) is exactly canceled by allowing for the 
slower “settling” times of very large particles. Equa¬ 
tions (4-7) of CW06, substituting Sc = 1 -|- St^ and 
Gett = ts -I- 1/St H from YL07, naturally lead to the same 
result: h/H = (at/St)^/^. While this has the proper 
behavior t o account for the parti cle layers in the sim¬ 
ulations of ICarballido et aP (1201 111 , which never exceed 
O.IH in thickness, it remains nonphysical in the small- 
St regime for reasons noted above, so must be converted 
using D95 Eq. (39) into a true particle scale height: 

hd/H= (l + St/at)-l/^ (B-2) 

which is valid for all combinations of St and at in 
global turbulence, where the energetic eddy frequencies 
are comparable to the orbit frequency. In nonturbulent 
nebulae with densely settled midplane layers, the near¬ 
midplane turbulence is shear driven and eddy times Te 
are faster, requiri ng modificat i on to these expressions; 
for discussions see ICuzzi et'H] (|1993l l or YL07 (unfortu¬ 
nately Fig. 6 of YL07 does not correctly represent the 
model of Cuzzi et al. 1993, which also includes a depen¬ 
dence on Te, leading to the ambigu ity of YL07 Fig. 7 
discussed bv ICarballido et ^ (1201111 . but this does not 
affect the overall results of YL07). 


APPENDIX C: CODE TESTS 
C.l. Coagulation 

Our numerical code uses the moments method of 
lEstrada fc Cuzzil (j2008ll to solve the Smoluchowski equa¬ 
tion ('Sec l2.4.1]) under the assumption that the dust dis¬ 
tribution can be treated as a powerlaw f{m) oc up 
to some fragmentation mass uif. The algorithm for the 
moments has been tested for both simple and realistic 
collisional kernels against the brute force solution to the 
coagulation equation for the cases in which the relative 
velocities between particles is systematic, or in the case 
where they are driven by turbulence. We specifically uti¬ 
lize the “modihed explicit” approach to follow the growth 
of the largest particle rp (Eq. (28), lEstrada fc Cuzzil 
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Figure 21. Upper panel: Comparison of the particle Stokes num¬ 
bers after 10^ years for four different choices for the number of bins 
per radius decade rip for the dust population. The blue curve is 
our nominal choice for this paper. Lower panel: As above, but for 
the migrator population. 

I2008D . This approach w as also compare d to the alterna¬ 
tive growth formalism of IGaraudI (1200711 and good agree¬ 
ment was found. Thus we do not repeat any tests of its 
accuracy here. 

C.2. Mass and Radius Binning 

In this section we check the influence of the choice of 
number rip of logarithmically spaced bins per decade ra¬ 
dius on the growth of particles in our code, both for the 
dust and migrator distributions. The moments method 
only requires a defined mass or radius histogram when 
calculating the collisional kernel, for example when solv¬ 
ing for the particle-to-particle relative velocities, which 
must be integrated over in order to calculate the growth 
rate of tol. This histogram is always well defined be¬ 
cause of the strict adherence to a powerlaw in mass or ra¬ 
dius. Nevertheless, because reasonably accurate growth 
requires sufficient sampling of the kernel quantities, es¬ 
pecially for the migrators, the number of bins chosen per 
decade can have an effect on the results. 

In the upper panel of Figure [m we compare the 
Stokes numbers Stu and St* (where defined) for a run 
of 10^ years for the fragmentation only case and an as¬ 
sumed at = 5 X 10“^, for four different choices of rip. 
The blue curve represents the value we have chosen for 
this paper. Even given the variation of Stu over orders of 
magnitude, all of the choices produce very similar results 
with the lowest value rip = 20 producing slightly more 
variation than higher values, but certainly within accept¬ 
able levels. The greatest difference occurs in the inner¬ 
most bin, which also happens to coincide with the edge 
of an EF. This is not surprising since the temperature 
variation is quite sensitive to the particle size distribu¬ 
tion. From these tests, it would appear that a choice of 
Up = 60 (whose curve appears most similar to the value 
of 100) may be a best choice going forward. It is also 



Figure 22. Temperature (black) and dust fraction (red) curves for 
different choices of the synchronization step. Our nominal choice 
is given by the solid curves. The relatively little variation indicates 
that the global processes such as the gas evolution and temperature 
calculation need not be done so frequently. 

possible to vary Tip with R. 

In the bottom panel of Fig. [211 we have instead varied 
the number of migrator radius bins per decade while us¬ 
ing our nominal value for the dust distribution. We com¬ 
pare our nominal value for the migrator bins rip = 100 
(blue curve) to smaller and larger values and find that 
this choice is more than sufficient as the rip = 100 and 
150 lie right on top of each other. We note that our 
treatment of the migrator distribution binning is different 
from that of the dust distribution which is effectively Eu- 
lerian as are all the usual solutions of the Smoluchowski 
equation, which lead to the well-known concerns about 
mass resolution (|Drazkowska et aTll2014l) . The migrator 
bins are, in fact, Lagrangian in that the migrator mass 
bins grow at different rates, and thus even a logarithmi¬ 
cally equally spaced grid would become unequal in its 
spacing with time. However, in order to maintain a con¬ 
sistency of spacing, we rebin the migrator distribution, 
strictly conserving mass, after every iteration. The gen¬ 
erally wide agreement between the variety of choices for 
rip give us confidence that growth is being treated well 
at all particle sizes we achieve in our models. 

C. 3. Variation of Ngync 

Global calculations in our code, such as the gas evo¬ 
lution, temperature and dust diffusion and advection do 
not need to be done every time step. This is the basis 
for our use of an asynchronous time stepping scheme. In 
this section we compare different choices of the number 
of steps fVsync our code executes before calling a syn¬ 
chronization step in which all properties of the nebula 
are brought to the same time. 

In Figure 1221 we compare four different choices of 
.^sync with our nominal choice of 100 shown by the solid 
curves. We plot both the temperature (black curves) and 
the dust fractional mass (red curves) for the same model 
as above after 10^ years. Again, we find that the agree- 
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ment across the different choices is more than adequate 
indicating that these processes are not changing quickly 
over a time scale of A^sync x Most of the variation 

is due to the sensitivity of the temperature to changes in 
the particle size distribution. Even these appear to be 
minor, especially in this case of high at in which growth 
can be rapid. 
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